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Abstract 



Serotonergic neurons of the dorsal raphe nucleus, with their extensive 
innervation of limbic and higher brain regions and interactions with the en- 
docrine system have important modulatory or regulatory effects on many 
cognitive, emotional and physiological processes. They have been strongly 
implicated in responses to stress and in the occurrence of major depressive 
disorder and other pyschiatric disorders. In order to quantify some of these 
effects, detailed mathematical models of the activity of such cells are required 
which describe their complex neurochemistry and neurophysiology. We con- 
sider here a single-compartment model of these neurons which is capable of 
describing many of the known features of spike generation, particularly the 
slow rhythmic pacemaking activity often observed in these cells in a variety 
of species. Included in the model are ten kinds of voltage dependent ion 
channels as well as calcium-dependent potassium current. Calcium dynamics 
includes buffering and pumping. In sections 3-9, each component is con- 
sidered in detail and parameters estimated from voltage clamp data where 
possible. In the next two sections simplified versions of some components 
are employed to explore the effects of various parameters on spiking, using a 
systematic approach, ending up with the following eleven components: a fast 
sodium current In a, a delayed rectifier potassium current Irdr, a transient 
potassium current I a, a low-threshold calcium current It, two high threshold 
calcium currents II and In, small and large conductance potassium currents 
Isk and Ibk, a hyperpolarization-activated cation current In, a leak current 
I Leak and intracellular calcium ion concentration Ca«. Attention is focused 
on the properties usually associated with these neurons, particularly long du- 
ration of action potential, pacemaker-like spiking and the ramp-like return 
to threshold after a spike. In some cases the membrane potential trajec- 
tories display doublets or have kinks or notches as have been reported in 
some experimental studies. The computed time courses of I a and It dur- 
ing the interspike interval support the generally held view of a competition 
between them in influencing the frequency of spiking. Spontaneous spiking 
could be obtained with small changes in a few parameters from their values 
with driven spiking. Spontaneous activity was facilitated by the presence of 
In which has been found in these neurons by some investigators. For rea- 
sonable sets of parameters spike frequencies between about 0.6 Hz and 1.2 
Hz are obtained. Anodal break phenomena predicted by the model are in 
agreement with experiment. There is a considerable discussion of in vitro 
versus in vivo firing behavior, with focus on the roles of noradrenergic input, 
corticotropin-releasing factor and orexinergic inputs. Location of cells within 
the nucleus is probably a major factor, along with the state of the animal. 
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Abbreviations 

5-HT, 5-hydroxytryptamine (serotonin); 5-HTP, 5-hydroxytrypto-phan; AHP, af- 
terhyperpolarization; acid; BK, big potassium channel; Ca;, internal calcium ion 
concentration; CB, calbindin-D28k; CBP, calcium binding protein; CDI, calcium- 
dependent inactivation; CNS, central nervous syytem; CR, calretinin; CRF, cor- 
ticotropin releasing factor; CSF, calcium source factor; D, duration (of spike); 
DA, dopamine; DRN, dorsal raphe nucleus; EPSP, excitatory post-synaptic poten- 
tial; FURA-2AM, Fura-2-acetoxymethyl ester; GABA, gamma-aminobutyric acid; 
HPA, hypothalamus-pituitary-adrenal cortex; HVA, high-voltage activated; ISI, in- 
terspike interval; LVA, low-voltage activated; mPFC, medial prefrontal cortex; PFC, 
prefrontal cortex; PV, parvalbumin; REM, rapid eye movement; SE, serotonin or 
serotonergic; SK, small potassium channel; SSRI, selective serotonin re-uptake in- 
hibitor; TEA, tetra-ethyl ammonium chloride; TPH, tryptophan hydroxylase; TTX, 
tetrodotoxin; VGCC, voltage-gated calcium channel. 
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1 Introduction: a summary of the properties of 
DRN SE neurons 

The last several decades have seen intensive experimental programs in neuroscience, 
endocrinology, psychiatry and psychology aimed at elucidating the involvement and 
responses of the many components of the nervous and endocrine systems in stress 
and to understanding their role in the complicated phenomenon of clinical depres- 
sion. Serotonergic neurons in the dorsal and other raphe nuclei, which extensively 
innervate most brain regions, have a large influence on many aspects of behav- 
ior, including sleep- wake cycles, mood and impulsivity (Liu et al, 2002; Miyazaki 
et al., 2011). With their influence on limbic and higher brain regions, they have 
important modulatory or regulatory effects on many cognitive, emotional and physi- 
ological processes (Lanfumey et al., 2005) including reward ( Nakamura et al., 2008; 
Kranz et al., 2010; Hayes and Greenshaw, 2011). They have been strongly impli- 
cated in responses to stress and in the pathophysiology of major depressive disorder 
and other stress-related psychiatric disorders (Neumeister et al., 2004; De Kloet et 
al., 2005; Firk and Markus, 2007; McEwen, 2007; Joels et al., 2007; Joels, 2008; 
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Linthorst and Reul, 2008; Savli et al., 2012). The therapeutical efficacy of selective 
serotonin re-uptake inhibitors indicates the key role of serotonin in depression. 

The involvement of serotonergic neurons arises through a variety of physiolog- 
ical, neurophysiological and endocrine processes. For example, the synthesis of 
serotonin, which may occur in the soma, synaptic terminals and to a lesser extent 
in some dendrites (Cooper et al., 2003; Adell et al., 2002) is from tryptophan which 
is supplied through the diet and enters the brain from blood. Within serotoner- 
gic cells the amino acid undergoes hydroxylation to 5-hydroxytrptophan (5-HTP) 
via tryptophan hydroxylase (TPH) which is synthesized only in the cell bodies of 
serotonergic cells. The activity of TPH is modifiable by some stressors and phar- 
macological agents. For example, through glucocorticoid receptors, sound stress 
increases the activity of TPH (Singh et al., 1990). 

Included in the targets of serotonergic neurons of the DRN are the hippocampus, 
amygdala, locus coeruleus, prefrontal cortex (PFC) and the paraventricular nucleus 
of the hypothalamus. Serotonergic input from the raphe into the hippocampus is 
important for regulating hippocampal neurogenesis (Balu and Lucki, 2009). 

Endocrine interactions and stress 

DRN SE neurons are known to be important mediators of endocrine (Chaouloff, 
2000; Amat et al., 2004) and behavioral (Maier and Watkins, 1998) responses to 
stressors. Bambico et al. (2009) found changes not only in spontaneous 5-HT neuron 
single-spike firing activity after chronic uncontrollable stress, but also changes in 
the number of spontaneously-active 5-HT neurons. Thus stress affects serotonergic 
neurotransmission by altering the firing rate of raphe 5-HT neurons, as well as 
the synthesis, release and metabolism of the transmitter and the levels of pre- and 
postsynaptic 5-HT receptors (Linthorst and Reul, 2005). 

The firing activity of DRN SE neurons is influential as serotonin (5-HT) is in- 
volved in the regulation of the hypothalamic-pituitary-adrenal (HPA) axis mainly 
by an effect exerted at a hypothalamic level (J0rgensen et al., 2002). The seroton- 
ergic activation of the HPA axis occurs by direct activation of CRF neurons in the 
paraventricular nucleus of the hypothalamus (Liposits et al., 1987). 

Stress may affect the serotonin system through corticotropin-releasing factor 
(CRF) fibers, which densely innervate the DRN where there are abundant CRF 
receptors (Pernar et al., 2004). DRN SE cells also have glucocorticoid receptors 
which are known to respond to high levels of Cortisol (Laaris et al., 1995). Kirby 
et al. (2000) demonstrated in vivo that a dense innervation of the DRN by CRF 
appeared to be topographically organized and showed that endogenous CRF in the 
DRN could alter the activity of 5-HT neurons. In vitro studies also suggest that 
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CRF modulates serotonergic neurotransmission (Lowry et al., 2000). The effects 
of CRF were principally excitatory and limited to a subpopulation of serotonergic 
neurons. Evidence that stress activates DRN 5-HT neuronal activity is indicated by 
an increase in 5-HT release within both the DRN and its projection areas (Lanfumey 
et al, 2005). 

Firing activity 

In most species, serotonergic neurons of the raphe nuclei have been found in vivo 
to have a slow, tonic pattern of firing. Usually spontaneous firing is found in rat 
DRN SE cells at rates around 0.5 to 2 Hz, as for example in Aghajanian et al. (1968), 
Aghajanian and Haigler (1974), Aghajanian and Vandermaelen (1982), Allers and 
Sharp (2003) with rates depending on the physiological state of the animal (Trulson 
and Jacobs, 1979; Jacobs and Azmitia, 1992; Urbain et al., 2006). There have 
been studies in freely moving rats in response to various stimuli (Waterhouse et al, 
2004). Some cells display bursting (Hajos et al., 1995, 1996; Schweimer et al., 2010). 
Aghajanian and Vandermaelen found that bursts of several spikes could be elicited 
with brief depolarizing currents of 0.5 to 1.0 nA. 

The robustness of the steady firing under a variety of conditions was taken to im- 
ply that spiking in serotonergic neurons was generated by intrinsic tonic pacemaker 
mechanisms. In fact, numerous intracellular recordings from dorsal raphe neurons 
show that spikes arise from gradual depolarizing ramps (Aghajanian and Vander- 
maelen, 1982; Vandermaelen and Aghajanian, 1983; Park, 1987) rather than the 
arrival of excitatory post synaptic potentials, although the latter could be present 
(Segal, 1985). 

DRN SE cells in vivo have a complex electrophysiological and neurochemical 
environment (Azmitia and Whitaker- Azmitia, 1995; Harsing, 2006; Lowry et al., 
2008) and it is not surprising therefore that recordings of spiking activity in vivo, 
usually made with chloral hydrate as anesthetic, have produced variable results. In 
slice, the activity of DRN SE neurons is usually reduced, which is usually attributed 
to the lack of excitatory synaptic input, which seems to contradict the claims that 
excitatory input is not required for their usual pacemaker-type firing. In Segal 
(1985), some cells required depolarizing currents of magnitude over 0.25 nA to fire. 
Interestingly, Vandermaelen and Aghajanian (1983) found that whereas in vivo 
most cells were spontaneously active (spiking), most cells were silent in slice - see 
also Kirby et al. (2003). Pan et al. (1990) recorded spontaneous spiking in 8 of 
104 neurons in slice. These findings contrast with earlier results by Mosko and 
Jacobs (1976) who claimed that the firing properties of raphe neurons in slice were 
mostly similar to in vivo, but the identification of the cells as serotonergic was not 
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attempted. With the application of either norepinephrine or the cti-adrenoceptor 
agonist phenylephrine, silent cells in slice may usually be made to fire regularly 
(VanderMaelen and Aghajanian, 1983; Judge and Gartside, 2006). In the former 
study, the injection of depolarizing currents of from 0.05 to 0.4 nA was also used to 
induce firing. 

Biophysical properties 

A survey of biophysical properties of rat DRN SE neurons reveals the following 
encapsulated view (for details see Tuckwell, 2012b). They have resting membrane 
potentials ranging from -75 to -58 mV with an average value (13 studies) of -64.4 
mV. Membrane time constants range from 7.4 to 53.5 ms, with a mean of 32.5 ms 
(8 studies) and input resistances from 30 to 600 MQ with a mean of 236 MQ (18 
studies). Such properties may vary across subfields (Calizo et al., 2011). Firing rates 
are variable, an average in vivo rate being 1.37 Hz (16 studies), usually with chloral 
hydrate, and ranging from 0.1 to 3.5 Hz. For 3 studies in slice, the average rate was 
1.07 Hz with range from 0.8 to 1.3 Hz. Spike durations are called long relative to 
many other types of neuron. It is difficult to give summary statistics for duration 
because of the many ways it is measured, but generally figures between 2 and 3.6 ms 
are cited, but sometimes up to 5 ms (Marinelli et al., 2004). Thresholds for firing 
are usually about 10 mV above resting potential and the afterhyperpolarization 
following spikes is on average 13.7 mV (9 studies). Spike amplitudes range from 
61 to 92 mV with a mean of 77.3 mV (9 studies), but again various definitions are 
employed. Estimated soma area, in /i 2 , is from 1073 to 2400 and estimated total 
somadendritic area is from 3914 to 11158 fi 2 . Capacitance can be estimated from 
area using the 1 /zF per square cm rule, and this gives a mean of 55 pF with a 
most likely value of 39 pF. Using instead the time constants and input resistances 
gives a higher mean of 89.5 pF. After considering all the data available on areas 
and capacitances, it was concluded that a typical DRN SE neuron has a soma 
dendritic area of 4000 fi 2 and a capacitance of 40 pF, these being the assumed 
values throughout most of the computations reported in this article - see Table 17. 
A very large cell might have an area of 9000 fi 2 and a capacitance of 90 pF. 

Inputs 

DRN SE neurons receive many inputs involving a variety of neurotransmitters 
(Jacobs and Azmitia, 1992) and the afferents are topographically organized (Lee 
et al., 2003; Hale and Lowry, 2011; Maximino, 2012), such as the important recip- 
rocal connections between the DRN and locus coeruleus (Kim et al., 2004). The 
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transmitters include norepinephrine and 5-HT itself. Norepinephrine, acting via 
ai-adrenoceptors, can accelerate the pacemaker activity of serotonergic neurons by 
closing potassium channels (Aghajanian and Sanders-Bush, 2002). 

5-HTi^ receptors (autoreceptors) in the 5-HT cell soma and proximal dendrites 
exert negative feedback control over DRN cell firing. Such serotonin-activated inhi- 
bition at autoreceptors may arise from volume transmission from self or neighboring 
DRN SE neurons as well as by recurrent and other axonal endings. A second pop- 
ulation of 5-HTia receptors is found throughout the dendritic field and modulates 
release of serotonin. Furthermore, 5-HTib/d receptors are found on the dendrites 
and preterminal axons where they modulate release of 5-HT (Stamford et al., 2000). 
Extrusion of 5-HT can occur through the reverse action of the serotonin transporter 
(Azmitia and Whitaker-Azmitia, 1995). 

Intracellular recordings from dorsal raphe neurons in brain slices show that 5- 
HTiyi agonists and 5-HT hyperpolarize the cell membrane, decreasing input re- 
sistance by opening K + channels, and decreasing high-threshold calcium currents 
(Penington and Kelly, 1990). A small component of this calcium current shows 
sensitivity to L-type calcium channel blockers and about 50% of what remains is 
sensitive to an N-type Ca 2+ channel blocker (Penington et al., 1991). There is 
evidence that the opening of K + channels via 5-HTi^ receptors in dorsal raphe neu- 
rons (Penington et al., 1993) is mediated by a pertussis-toxin-sensitive G protein 
(Williams et al., 1988). There are also glutamatergic inputs directly on DRN SE 
neurons, from the PFC, lateral habenula, the hypothalamus and other regions, as 
revealed by glutamate transporter tracing (Soiza-Reilly and Commons, 2011). 

The mPFC not only receives strong inputs from serotonergic neurons in the 
DRN but also sends projections to this nucleus. Electrophysiological investigations 
in the rat DRN reveal that most serotonergic neurons are inhibited by electrical 
stimulation of the mPFC, suggesting that this pathway is more likely to synapse 
onto neighboring GABA-ergic neurons rather than onto 5-HT cells (Jankowski and 
Sesack, 04). The mPFC thereby exerts control over the DRN and that may be in- 
volved in the actions of pharmaceutical drugs and drugs of abuse (Goncalves, 2009). 
Thus, in addition to feedback systems involving separate nuclei and structures, local 
networks are found within the dorsal raphe nucleus involving interactions between 
5-HT and local inhibitory GABA-ergic and excitatory glutamatergic neurons (Agha- 
janian and Sanders-Bush, 2002) 

Perspective 

Mathematical modeling of neurophysiological dynamics has been pursued for 
many different nerve and muscle cell types. Some well known neuronal examples 
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are thalamic relay cells (Huguenard and McCormick, 1992; Destexhe et al., 1998; 
Rhodes and Llinas, 2005), dopaminergic cells (Komendantov et al., 2004; Putzier 
et al., 2009; Kuznetsova et al., 2010), hippocampal pyramidal cells (Traub et al., 
1991; Migliore et al., 1995; Poirazi et al., 2003; Xu and Clancy, 2008), neocortical 
pyramidal cells (Destexhe et al., 2001; Traub et al, 2003; Yu et al., 2008) and 
molluscan interneurons (Vavoulis et al., 2010). Cardiac myocytes have also been 
the subject of numerous computational modeling studies with similar structure and 
equivalent complexity to that of neurons (Faber et al., 2007; Williams et al., 2010). 

It is a complicated task to model all the details of an SE neuron. Some authors 
have focused on models of synaptic release (Best et al., 2010) including the effects 
of SSRIs (Geldof et al., 2008). However, the dynamics of the various currents are 
important in determining the sequence of action potentials which will determine, 
despite the many complexities of the axonal branchings, the amounts of transmitter 
release. Hence the focus here is on the modeling of the process of spike generation 
which is fundamental as a prelude to understanding the responses of these cells 
to activation of both neurophysiological (including synaptic) and endocrinological 
receptors, as well as the effects of this spiking and subsequent 5-HT release at 
numerous sites. 

The firing patterns of these cells have been much studied and many properties of 
the ionic currents underlying action potentials have been investigated (Aghajanian, 
1985; Segal, 1985; Burlhis and Aghajanian, 1987; Penington et al., 1991, 1992; 
Penington and Fox, 1995, Chen et al., 2002). In the present work, modeling of spike 
generation in DRN SE neurons has emphasis on intrinsic properties, being a first 
step in the construction of a quantitative theory of how these neurons are influenced 
by their many inputs and how they influence their numerous target cells. 

According to Aghajanian and Sanders-Bush (2002) and several others, the pace- 
maker activity of serotonergic neurons results from the interplay of several intrinsic 
ionic currents including a voltage-dependent transient outward potassium current 
I a, a low-threshold inward calcium current It, and a calcium-activated outward 
potassium current. There have been conflicting ideas about the role of Ia (Burlhis 
and Aghajanian, 1987; Segal, 1985). Since the basis of pacemaker-like activity in 
these cells is not yet fully understood, one of the aims of the present article is to 
investigate these purported mechanisms in a quantitative way. 
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2 Membrane currents 



For a single compartmental neuron model, the differential equation for the mem- 
brane potential is, 

C^ = -I,V(0) = V , (1) 

where / is the sum of all current types and Vq is the initial value of V, taken in 
most of what follows to be the resting membrane potential, Vr. Here depolarizing 
currents are negative, V is in mV, and t is in ms so that if / is in nA, then the 
membrane capacitance is in nF. If the i-th component of the current is denoted by 
Ii so that / = ^2 Ii, then following the Ho dgkin- Huxley (1952) model, generically 
and generally, each component current, ij, is taken as a product of activation and 
inactivation variables, a maximal conductance, gi t max, and a driving force which is 
V — Vi where Vi is usually at or near the Nernst equilibrium potential. 

For noninactivating currents there is an activation variable m raised to a certain 
power p > 1, not necessarily an integer, so that 

I i = g i , ma xrn?(V-V i ). (2) 

If the current inactivates, then the current contains an inactivation variable h which 
is usually raised to the power 1 so 

k = 9i, ma xm p h(V - Vi). (3) 

Sometimes the activation or inactivation variables depend on, or also on, calcium 
ion concentration. A constant field expression may be used instead of the linear 
term V — Vi - see Section 4 on calcium currents. 

Activation and inactivation variables are determined by differential equations 
which are conveniently written in the forms 

dm moo — rn , , s 
dt r m 

dh _ hoc - h , s 

dt ~ r h [ ' 

where and hoo are steady state values which depend on voltage. The quantities 
r m and Th are time constants which may also depend on voltage and/or calcium 
concentration. We adhere to the convention throughout that all parameters as 
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written are positive. For accounts of the dynamics of many types of current, see 
Levitan and Kaczmarek (1987) and Destexhe and Sejnowski (2001). 

Table 1 gives the notation for the 11 membrane currents included in the present 
model, so that 

dV 

C— = -[I A + IkDR + I M + I T + h + In + I H + iNa + I B K + hK + heak], (6) 

though in the computations various components are sometimes omitted. In addition 
there is a differential equation describing the evolution of the internal calcium ion 
concentration in terms of sources and a pump, 

dCoj 

r, J sources J pump- 



Table 1: Notation for membrane currents 



Symbol 



Description 



Voltage-gated potassium 

I A 

Irdr 
Im 

Calcium 

It 
II 
In 

Calcium-acivated potassium 

Ibk 

ISK 

Hyperpolarization activated cation current 

Ih 

Fast transient sodium current 

I Na 

Leak current 

I Leak 



Fast transient A-type 
Delayed rectifier potassium 
M-type potassium 

Low threshold calcium, T-type 
High threshold calcium, L-type 
High threshold calcium, N-type 

Large conductance channels 
Small conductance channels 
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3 Voltage-dependent potassium currents 



There is evidence for three types of voltage-gated potassium currents in DRN SE 
neurons. These are the transient I a, and two which are usually considered to be 
non-inactivating or very slowly inactivating, being the delayed rectifier, Irdr, and 
the M-current, Im- 

M-type potassium current, Im 

The M-current is associated with the potassium channel types K v 7.1 to K v 7.5 or 
the KCNQ family (Gutman et al., 2005). Such a current activates slowly and inacti- 
vates extremely slowly, more so than the usual delayed rectifier. M-type potassium 
currents are considered to play a significant role in adaptation (Storm, 1990; Benda 
and Herz, 2003). M-type potassium currents have been included in several com- 
putational models of neurons, especially hippocampal and cortical pyramidal cells 
(Lytton and Sejnowski, 1991; Poirazi et al., 2003; Vervaeke et al., 2006; Xu and 
Clancy, 2008). However, such currents are associated with extremely small conduc- 
tance densities, being of the order 1000 to 4000 times less than that of the usual 
delayed rectifier. 

The potassium M-current is written as 

Im = gM,maxm M (V ~ V M ) (7) 

where Vm is the reversal potential, taken as around the Nernst potential, Vk, for 
potassium. The activation variable m M satisfies an equation like (4) and the steady 
state is written 

1 

mM '°° _ i + e -(v-v Ml )/k Ml ( 8 ) 

with a time constant 

T m M = T m,Mci (9) 

a constant. In voltage clamp experiments on potassium currents similar to those 
in our recent report on I a (Penington and Tuckwell, 2012), in some cells we have 
observed a persistent current which seems to have properties similar to an M-current. 
An example is shown in Figure 1 where the current is decomposed into a component 
likely to be I a, which declines exponentially, and a persistent component which is 
putatively an M-current with time constant of activation of about 50 ms and no 
sign of inactivation. 



12 



1600 



< 



1400 - 



1200 - 



1000 - 



200 - 



800 - 



600 - 



400 - 



- 




-200 



50 



100 



150 



200 



250 



300 



Time (ms) 



Figure 1: Decomposition of current under voltage clamp (-60 mV to 50 mV) in 
a putative DRN SE neuron. The exponentially decaying curve is likely to be I a 
whereas the persistent current is possibly Im- 

Parameter values for Im 

Half- activation potentials for Im have been obtained in subtypes expressed in Chi- 
nese hamster ovary cells and ranged from -28.7 mV to -11.6 mV with corresponding 
slope factors from 10.1 to 12.9 (Tatulian et al., 2001). In various modeling and 
other studies the half- activation potentials have ranged from as low as -50 mV to 
-20 mV with slope factors around 9 or 10. Time constants for activation have been 
from about 50 ms to 200 ms, depending in some cases on V. The following set of 
parameters (Table 2) may be chosen for Im, where Vk is the Nernst potential for 
potassium. 



Transient potassium current, I a 

The transient potassium current I a was documented in early studies of DRN SE 
neurons by, for example, Aghajanian (1985), Segal (1985) and Burlhis and Agha- 



Table 2: Standard set of parameters for Im 



Parameter Value Parameter Value 
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janian (1987). I\ has been posited to play an important role in determining the 
cell's firing rate. However, there have been conflicting ideas about the role of I a 
(Burlhis and Aghajanian, 1987; Segal, 1985). Six kinds of K v channels have been 
linked to Ia (Gutman et al., 2005) but according to Serodio and Rudy (1998) and 
other studies, in the mammalian brain, somatodendritic Ia currents are likely to be 
carried by channels of the Shal family K v 4.1, K v 4.2 and K v 4.3. Using in situ hy- 
bridization histochemistry, Serodio and Rudy (1998) found that in the DRN, K v 4.3 
signals were strong, K v 4.2 weak and K v 4.1 below background. However, Pearson 
et al. (2006) in their study of DRN and locus coeruleus, found evidence for K v 4.2 
in the DRN of adult rats. In a recent report (Penington and Tuckwell, 2012), we 
have examined various data which suggested that the channels carrying I a in the 
putative serotonergic neuron analyzed are of the K v 4.2 type, but there was insuf- 
ficient evidence to completely rule out the other two Shal members. In support 
of this idea are the findings of Bekkers (2000) who postulated that in rat layer 5 
cortical pyramidal cells, I a was carried by K v 4.2 channels, and of Kim et al. (2005), 
who found that K v 4.2 is the main contributor to Ia in hippocampal CAl pyramidal 
cells, determining several basic characteristics of the action potential, including its 
width. 

In their model of thalamic relay cells, Huguenard and McCormick (1992), em- 
ployed two types of Ia with somewhat different kinetic properties. However, we do 
not have evidence for any but one type which we denote simply by I a with the form 

Ia = gA,maxm A h A (V - V A ) (10) 

where tua and h A satisfy equations like (6) and (7). The power 4 is used for itla 
based on our previous observations (Penington and Tuckwell, 2012), whereas other 
authors have used 3. The steady state activation is written 

1 

mA <°°- i + e -(v-v Al )/k Al ( n ) 

with a time constant 

TmA =aA + cosh((V - A V A2 )/kA 2 y (12) 

where there are 4 parameters a a, Va 2 , kA 2 - The steady state inactivation variable 
is set to 

hA '°° = i + e (v-v A[i )/k A3 ( 13 ) 
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with corresponding voltage-dependent time constant, 



T m A — C A "+ 



d A 



(14) 



cosh((v-v A4 )/k A4 y 



Parameter values for I A 

The standard set of kinetic parameters for I A is based on one given in Huguenard 
and McCormick (1992), with suitable approximations for the time constants. In the 
case of the time constant for inactivation, the approximation is valid over values of 
V which are likely to occur during the activity of a DRN SE neuron. The use of 
these parameters is validated by the results obtained below in simulation of some 
voltage clamp experiments of Burlhis and Aghajanian (1987) on DRN SE neurons 
in slice. 



Delayed rectifier potassium current, Irdr 

Evidence for a delayed rectifier potassium current in DRN SE neurons has been given 
by Penington et al. (1992) and Liu et al. (2002) who found it was substantial in 
the repolarization phase of action potentials. There are apparently no explicit data 
on the properties of Irdr i n these cells, but since the original modeling of squid 
axon action potential by Hodgkin and Huxley (1952) there have been numerous 
approaches to a quantitative framework. 
We take the form for this current to be as 



Table 3: Standard set of parameters for I A 



Parameter Value Parameter Value 




IkDR — 9KDR,maxn 



nk {V-V KDR ) 



(15) 
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where n (the traditional symbol) is the activation variable which satisfies a differ- 
ential equation like (4), Vkdr is the reversal potential, and n& is usually an integer 
between 1 and 4, inclusively. The steady state activation is written 

n °° = i _|_ e -(v-v KDRl )/k KDRl ( 16 ) 



and the time constant 



okdr m „s 
&KDR H TT7T} 77 wl 7 I 1 'J 

cosh((V - V KDR2 )/k KDR2 ) 



Parameter values for Ikdr 

Table 4 gives parameters on the activation used in various studies (some of which 
are in conjunction with experiment) and it can be seen that all integer powers from 
1 to 4 have been used for n, with 4 the value in the squid axon model (Hodgkin and 
Huxley, 1952). We will choose the power 1 in the absence of specific data as the 
properties are then simpler to relate to the steady state activation function, but also 
examine other values with appropriate changes in other parameters. It can be seen 
that there is considerable variability assumed or measured in the kinetic properties 
of Ikdr in various neurons. Note that Traub et al. (1991) values were also used by 
Migliore (1995) and Xu (2008) and that in Bekkers (2000) the index may vary with 
V. 



Table 4: Examples of parameters for activation of Ikdr 



Source 


Form 


Vkdr 1 


kKDRi 


V Rest 


Train T rnax 


V at T max 


Traub (1991)* 


n 


-18.2 


10 


-60 


0.8-4.2 


-31 


Belluzzi (1991) 


n 


-6.1 


8 


-75 


2.5-25 


-27 


Schild (1993) 


n 2 


-7.2 


11.8 


-51 


4-64 


-32 


Bekkers (2000) 


TV 


-9.6 


13.2 




1.5-20 


-20 


Traub (2003) 


4 

n 


-29.5 


10 


-62 to -68 


0.25-4.6 


-10 


Poirazi (2003) 


n 2 


-42 


2 


-70 


2.3-3.5 




Molineux (2005) 


n 


-35 


4 








Komendantov (2007) 


n 3 


-18.3 


10 


-58 


1-2-6.3 


-34 


Anderson (2010) 


n 


-25 


4 


-65 to -60 
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Based on the properties of the repolarization and afterhyperpolarizations in 
spikes in DRN SE neurons, it is postulated that the activation time constant for 
Ikdr would not be very large (less than 5 ms) and that it would be mainly activated 
between about -40 mV and + 20 mV. Thus, considering the various values given 
in Table 4, the following parameters (Table 5) are taken as the standard set. The 
reversal potential is taken to be Vk- 



Table 5: Standard set of parameters for Ikdr 



Parameter 


Value 


Parameter 


Value 


Vrdr 1 


-15 


bKDR 


4 


k-KDR 1 


7 


Vrdr 2 


-20 


&KDR 


1 


kxDR 2 


7 



To estimate the conductance associated with the delayed rectifier potassium cur- 
rent we appeal to some voltage clamp activation experiments done with and without 
TEA on dissociated cells. The methods were described in Penington and Tuckwell 
(2012). Results for one DRN serotonergic cell are shown in Figure 2. From the 
uppermost recordings (clamps from -110 mV to mV), the final current amplitude 
without TEA minus the final amplitude with TEA gives 711 pA. Mathematical 
modeling of the potassium delayed rectifier current then yields an associated con- 
ductance of 0.0085 fiS and since a dissociated cell has a capacitance of about 20 pF 
(Penington and Fox, 1995), this translates to 0.425 nS/pF, at room temperature. 
Allowing for an increase due to an increase in temperature, a figure of 0.641 nS per 
pF seems reasonable for the maximal conductance of Ikdr at body temperature. 

4 Calcium currents 

Calcium currents, which are found in all excitable cells, have been generally di- 
vided into the two main groups of low-threshold or low- voltage activated (LVA) and 
high-threshold or high- voltage activated (HVA). The former group contains only 
the T-type (T for transient) and the latter group consists of the types L, N, P/Q 
and R (L for so called long-lasting, N, either for neither T nor L, or neuronal, P 
for Purkinje, and R for resistant). Calcium channels have up to four subunits, a±, 
Oi2-5, (5 and 7, which may exist in different forms, and which modulate the con- 
ductance and dynamical properties of the channel (Dolphin, 2006, 2009; Davies et 
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Figure 2: Current versus time with voltage clamp in activation experiments with 
TTX alone (left column) and TTX + 20 mM TEA (right column). Steps from -120 
mV to a maximum of mV. Markers 100 ms and 100 pA. 

al., 2010). The ten forms of the conducting pore subunit, a±, lead to an expansion 
of the above groups (Catterall et al., 2005; Dolphin, 2009) to 10 main subtypes. 
According to the accepted nomenclature, L-type channels consist of the subtypes 
Ca v l.l — Ca v 1.4. The remaining "high-threshold" currents, P/Q, N and R are re- 
spectively Ca v 2.1 -Ca v 2.3 and the T-current subtypes are Ca v 3.1 — Ca v 3.3. Within 
the subtypes, various configurations of other subunits lead to channels, with quite 
different properties (Dolphin, 2009). One cannot therefore ascribe definite parame- 
ters in the dynamical description of calcium currents based on the subtypes. L-type 
channels, for example, display a wide variety of activation properties (Tuckwell, 
2012a). The properties of channel types may not be the same in all locations of the 
same cell. For example, somatic N-type channels do not inactivate in rat supraoptic 
neurons but do in synaptic terminals of the same cell and in other cell types (Joux 
et al., 2001). 

Calcium T-type current, It 

The T-type (low threshold) calcium current has been implicated in pacemaker ac- 
tivity in some cells and also in bursting activity (Perez- Reyes, 2003; Catterall et al., 
2005) and amplification of dendritic inputs (Crandall et al., 2010). The three main 
subtypes are designated Ca v 3.1, Ca v 3.2 and Ca v 3.3 or aiG, a\H and a±I, respec- 
tively. Since the pioneering experimental studies of Segal (1985) and Burlhis and 
Aghajanian (1987), spike generation and the pacemaker cycle in DRN SE neurons 
have been posited to be triggered by T-type currents. See also Jacobs and Azmitia 
(1992) and Aghajanian and Sanders-Bush (2002). 
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Mathematical expressions for I T 

The mathematical form for the T-type current has usually been chosen to be the 
linear expression (Destexhe et al., 1994, 1996; Rybak et al., 1997; Amini et al., 1999; 
Sanchez et al., 2003), 

l T = g T ^ max m^h T (y — V C a t rev), (18) 

although some authors have m T rather than m T (Rhodes and LLinas, 2005). 
and Jit satisfy differential equations like (4) and (5). The advantage of this form for 
It is that the associated conductance can be compared with other ion channels for 
which the linear form is usually employed - see the simplified model in section 10. 
However, because the linear form may sometimes lead to inaccuracies especially for 
calcium ion currrents and in particular for low internal calcium ion concentrations 
(Belluzzi and Sacchi, 1991; Huguenard and McCormick, 1992; Poirazi et al., 2003), 
it is often desirable to use the constant field form 



IV - 



TCaj — Ca G e v <> 1 
I T = k T m T h T V± — — j - (19) 

1 - e v o 

where 

4F 

k T = 1000Ap T P T —. (20) 

Here A is membrane area in sq cm, is the channel density, P T is the single 
channel permeability in cm 3 per second, Ca;, Ca Q are the (time dependent) inter- 
nal and external concentrations of Ca 2+ in mM, F is Faraday's constant (96000 
Coulombs/Mole), 

Vo = ^r, (21) 

and the factor 1000 converts the current to nA. The time constants for activation 
and inactivation are denoted by r mT and Th T . 

For current through some Ca 2+ channels, there is not only volt age- dependent 
inactivation but also inactivation by internal calcium, called calcium dependent 
inactivation (CDI). When present, CDI can exert a considerable effect (Tuckwell, 
2012a). However, for the low threshold current It, there is no CDI (Budde et al., 
2002; Dunlap, 2007), but it was included in a model of a CAl pyramidal cell (Poirazi 
et al., 2003). We will assume there is no CDI for I T in DRN SE neurons. 
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Voltage clamp data: m Too , h T ,oo 

Voltage clamp data are available for activation and inactivation of It in various cell 
types and some data are reviewed in Perez- Reyes (2003) and Catterall et al. (2005). 
For the steady state activation we put 

mT '°° = i + e -(v-v Tl )/k Tl > ( 22 ) 

where Vt x is the half- activation voltage and and fc^ is the slope factor. Similarly 
for the steady state inactivation 

hT '°° = i + e (v-v T3 )/k T3 ■ (23) 

Table 6 gives values of the constants V^, hr l , Vr 3 and /ct 3 obtained for some 
CNS cells. All experimental data are for dissociated cells. There are also data 
(not shown) for the three subtypes a\G, a\H and a\I , expressed in HEK-293 cells 
(Klockner et al., 1999). The half-activation potentials ranged from -45.8 mV to 
-43.8 mV, the half-inactivation potentials ranged from -72.8 mV to -72.0 mV and 
the slope factors for inactivation ranged from -4.6 mV to -3.7 mV. In their modeling 
of thalamic relay cells, Rhodes and Llinas (2005) used the data of Huguenard and 
McCormick (1992), shifted in the direction of depolarization by 5 mV to prevent a 
large It at rest. Destexhe et al. (1998) also shifted the same data by a few mV in 
the hyperpolarizing direction to allow for different external Ca 2+ concentrations and 
in the depolarizing direction as well to obtain agreement with their experimental 
observations. 



Table 6: Parameters for steady state activation and inactivation for I T 



Neuron type 


v Tl 


k Tl 


V T3 


k n 


Thalamic relay, Huguenard (1992), 


-57 


6.2 


-81 


-4 


Thalamic relay, Huguenard & Prince (1992) 


-59 


5.2 


-81 


-4.4 


Thalamic reticular, Huguenard & Prince (1992) 


-50 


7.4 


-78 


-5.0 


Thalamic relay, Destexhe (1998) 


-56 


6.2 


-80 


-4 


Thalamic relay model, Rhodes (2005) 


-63 


6.2 


-80 


-4 
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Table 7: Maximal currents for It in a DRN SEneuron : V from -80 to V\ 



v 1 


Ifnax 


I-max 


Ifnax 


(mV) 


Expt (pA) 


V Tl = -57ml/ 


V Tl = -44ml/ 


-35 


230 


230 


230 


-40 


161 


228 


160 


-45 


93 


208 


83 


-50 


44 


165 


32 



For presumed DRN SE neurons (dissociated), some voltage clamp data for I T 
were given by Penington et al. (1991). For a holding potential of -80 mV, records 
of It were given for test potentials of V\ = —35, —40, —45, —50 mV and estimates 
of the peak currents are given in column 2 of Table 7. Using the kinetic data of 
Huguenard and McCormick (1992), with V Tl = —57 mV, gave the peak currents in 
column 3 which are not close to the experimental values. Only if V Tl was shifted to 
-44 mV, with or without a similar shift in the inactivation curve, was it possible to 
obtain an approximate agreement with the experimental data. Such discrepancies 
could be due to a number of factors such as temperature, differing compositions 
of intracellular and extracellular fluids or that in the experiments 5 mM Ba 2+ was 
the charge carrier. However, many reports (Huguenard and Prince, 1992; Bourinet 
et al., 1994; Rodriguez-Contreras and Yamoah, 2003; Durante et al., 2004; Goo et 
al., 2006) indicate that barium shifts the kinetics in the hyperpolarizing direction 
relative to calcium, although in one study a depolarizing shift was noted (Lipscombe 
et al., 2004). 



Time constants 

The corresponding time constants have been written in several ways including a 
constant plus the reciprocal of the sum of two exponentials (McCormick and Hugue- 
nard, 1992; Destexhe et al, 1994; Rybak et al. 1997; Rhodes and Llinas, 2005) and 
a Gaussian (Amini et al., 1999). For activation the following form was able to fit 
well several of the functions employed 

T mT =a T + — bT (24) 
cosh((V - V T2 )/k T2 ) 

where there are 4 parameters ax, br, Vr 2 , kx 2 - 

For the inactivation time constant Th T the same forms were used as for r mT by 
Destexhe et al. (1994), Rhodes and Llinas (2005) and Amini et al. (1999), and 
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a function with discontinuous derivative was used by McCormick and Huguenard, 
(1992) and Rybak et al. (1997). A good fit was found possible with the form used 
by Amini et al. (1999), namely 

r hT = Ct + d T e-( (v - VT ^ kT ^ 2 , (25) 

where there are 4 parameters ct, dx, Vr 4 , k Ti . In some calculations, such as those 
mentioned in the previous subsection, the forms given by Huguenard and Mc- 
Cormick (1992) were employed. 

Voltage clamp results of Burlhis and Aghajanian (1987) 

Where possible, in order to obtain estimates of the properties of various component 
currents, experimental data can be compared with theoretical predictions. Burlhis 
and Aghajanian (1987) performed several experiments in brain slice, designed to 
examine quantitatively the currents underlying pacemaking in intact rat DRN SE 
neurons In one such experiment, called an "anodal break" (see also Section 12) , a 
single-electrode voltage clamp was applied with a depolarizing step from -80 mV to 
-56 mV as seen in their Figure 2B. 
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Figure 3: Steady state activation and inactivation functions versus membrane po- 
tential for I a and It with parameters from Huguenard and McCormick (1992). Here 
are plotted m Too and m\ OQ . 
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With the start of the depolarizing step, a transient outward current was observed 
with a peak magnitude of about 615 pA at around 7 ms. This was followed by a 
slower but transient inward current of peak magnitude about 92 pA around 67 ms. 
The early component was posited to be the fast transient potassium current Ia 
and the inward current was the transient It calcium current, being the component 
responsible for what Burlhis and Aghajanian (1987) called the prepotential which 
occurred before an action potential. These two component currents were modeled 
under the voltage clamp, -80 mV — > -56 mV, using the parameters of activation 
and inactivation, including steady state activation and inactivation functions and 
time constants given for thalamic relay cells given by Huguenard and McCormick 
(1992). Note that the latter cells have resting potentials very similar (about -63 
mV) to those DRN SE neurons. 
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Figure 4: Model results for I a + It fo r comparison with the voltage clamp results, 
-80 mV — > -56 mV, of Figure 2B in Burlhis and Aghajanian (1987). 

In this modeling, the constant field form for It was used. The steady state values 
of m T oo and m\ and also /i-r,oo and /ia,oo are shown in Figure 3 where it can be 
seen that the half-activation values are about -56 mV for It and about -48 mV for 
I a- A computed result is shown in Figure 4. The top part shows the voltage step 
and the lower part shows the time course of the sum of I a and It- The maximal 
conductances had to be estimated by trial and error which gave approximately for 
Ia, the value gA,max = 0.479 /iS. The estimated value of /c^ defined by Equ. (20) 
which resulted in the Ia + It curve in Figure 4 was fc^ = 0.0825. 
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Table 8: Conductance densities for I a for various cells 



Source 


Cell type 


9A,max 


(1st auth) 




( /iS cm" 2 ) 


Lytton (1991) 


Cortical pyramid 


100000 


Amini (1999) 


Midbrain dopaminergic 


445 


Athanasiades (2000) 


Medullary respiratory 


6730 


Traub (2003) 


Cortical pyramid 


30000 


Komendantov (2007) 


Magnocellular neuroendocrine 


15000 to 20000 


Saarinen (2008) 


Cerebellar granule 


445 


Penington (2012) 


DRN 


1600 


This study 


SE DRN 


12000 



The estimate of gA,max translates to a density of about 12000 /xS per cm 2 , on the 
assumption of a soma-dendritic area of 4000 /im 2 . This figure is compared with 
those from some other neurons in Table 8 where an extremely broad range of values 
is seen, from a smallest value of 445 /xS per cm 2 to 100000 /xS per cm 2 . With 
the estimated value gA,max from the voltage-clamp measurements of Burlhis and 
Aghajanian (1987), the predicted time course for Ia + I? shown in Figure 4 is in 
broad agreement with the experimental results. However, this is only a heuristic 
result which may be useful as a quantitative guide, because of the unknown effects 
of various other membrane currents and fluxes which might affect the experimental 
results. 

In a second voltage clamp experiment, Burlhis and Aghajanian (1987) used 
CsCl-filled electrodes to suppress I a- The holding potential was again -80 mV 
and the peak inward current was plotted against test potential in their Figure 5D. 
According to Puil and Werman (1981), Cs + blocks not only 1a but also Irdr and 
calcium-activated potassium currents (but see also Sanchez et al., 1998). A model 
of an DRN SE neuron with all potassium currents blocked was employed to examine 
the current- voltage relation and the results are shown in Figure 5. There is broad 
agreement with the experimental result as there is a fairly sharp upswing in the 
peak inward current from about -70 mV to -50 mV, indicating that It has been 
switched on at about the right voltages to give the prepotentials which trigger spikes 
in these cells. 
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Figure 5: An attempt to emulate voltage-clamp results of Burlhis and Aghajanian 
(1987, Figure 5D). From a holding potential of -80 mV, test potentials gave rise to 
computed peak inward currents of magnitudes shown, indicating broad agreement 
with the experimental results. 

Parameter values for It 

The standard set of kinetic parameters for It given in Table 9 is mainly taken 
from Huguenard and McCormick (1992), with suitable approximations for the time 
constants. In the case of the time constant for inactivation, the approximation is 
valid over values of V which are likely to occur during the activity of a DRN SE 
neuron. 



Table 9: Standard set of parameters for It 



Parameter 


Value 


Parameter 


Value 


v Tl 


-57 


V T3 


-81 m 


k Tl 


6.2 


k Ta 


-4 


OjT 


0.7 


c T 


28 


b T 


13.5 


cIt 


300 


V T2 


-76 


v n 


-81 


k T2 


18 




12 
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Calcium L-type current, II 

There are four main subtypes of L-type calcium currents, designated, according to 
the accepted nomenclature, Ca v l.l — Ca v 1.4. In most central nervous system cells 
the subtypes are mainly Ca v 1.2 or Ca v 1.3. According to a recent review (Tuckwell, 
2012a), L-type calcium currents often inactivate by both voltage-dependent and 
calcium dependent mechanisms (CDI). The general form for this current in the 
linear version is similar to that of It in Equation (18) with L replacing T throughout. 
Similarly with the constant field formalism the current is 

r -2Zn 

9 Ca ; - Ca Q e 
I L = k L m 2 L h L f L V. ^ 4^ (26) 

1 - e v o 

where most of the symbols are as defined for I T except that now f L represents the 
CDI. rriL and satisfy differential equations like (4) and (5). For neurons the 
powers 1 and 2 are about equally frequently employed for m^. The equation for fi 
is 

* - r f ' (27) 

where /z, i00 (Cai) is the steady state value /z,(Caj, oo) and 77 is a time constant which 
may depend on Ca;. Here /l iOG is defined as 

/l,oo — - ^ ^Ca;^" (^) 

where Kf, in mM, is the value of Ca; at which the steady state inactivation has half- 
maximal value. Knowledge of 77 is scant, and one can adopt the usual assumption 
that it is very small (especially relative to t^l) so that fi is set equal to its steady 
state value. We will also adopt the usual value n — 1 (Standen and Stanfield, 1982). 
Thus we may put 
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, 1 [Cai - Ca Q e y o 1 
h = k L m\h L V. L-L (29) 

1 + \-kJ) 1 - e v 

Assumptions about the parameters for II are not expected to have grave conse- 
quences for the spiking dynamics because the contributions from the L-type current 
are presumed to be less than about 10% of the whole cell calcium currents, based on 
the observations of Penington et al. (1991) that the L-type contribution was small 
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(about 4%) in dissociated cells and that the threshold for activation is not likely to 
occur around potentials near the threshold for action potentials. 
For the steady state L-type activation we put 

m ^ = i + e -(v-v Ll )/* £l ( 3 °) 

and 

T r n L =aL-\ TJ7T} , \r x T, 7) (31) 

cosh((y + V L2 )/k L2 ) 

where there are 4 parameters a^, &l, Vl 2 , fcj, 2 . These two functions, m i oo and r„ tL 
were well-fitted, in accordance with experimental data, to the more complicated 
forms in Athanasiades et al. (2000) and Komendantov et al. (2007). The same 
remarks apply to the steady state voltage-dependent inactivation variable which is 
set to j 

hL >°° = i + e (v-v La )/k La ( 32 ) 
with a voltage- independent time constant, 

7~h,L = const. (33) 



Parameter values for I L 

For DRN SE neurons there are no explicit voltage clamp data from which kinetic 
parameters for activation and inactivation can be estimated for I L . In Penington et 
al. (1991), single channel currents (slope 23 pS) were given for L-type with voltages 
ranging from about -20 mV to + 30 mV. It is likely that the L-type current in these 
cells belongs to the high neuron group summarized in Table 3 of Tuckwell (2012a). 
Aghajanian (personal communication) found that L-channel blockers had little effect 
on firing rates, which implies that these channels are not of the pacemaking type as 
they are in dopaminergic neurons (see Discussion), in which the activation occurs 
at much lower voltages, as in the low neuron group of Tuckwell (2012a). 

The average data is employed as the standard set in Table 10 in the present 
article. 



Calcium N-type current,/^ 

In voltage clamp experiments on DRN SE neurons, Penington et al. (1991) found 
evidence that among high-threshold calcium currents, N-type channels contributed 
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Table 10: Standard set of parameters for II 



Parameter 


Value 


Parameter 


Value 


v Ll 


-18.3 


ki 2 


15 




8.4 


Vl 3 


-42 


a L 


0.5 




13.8 


b L 


1.5 


r h L 


200 


Vl 2 


-20 







about 40% of total calcium current in dissociated cells. Most of the remaining frac- 
tion was not blocked by cv-conotoxin and was called N-like. These two components 
are here lumped together as N-type. Whether CDI is appreciable for N-type cal- 
cium currents is uncertain (Budde et al., 2002), although some evidence for it has 
been obtained in chick dorsal root ganglion neurons (Cox and Dunlap, 1994) and 
rat sympathetic neurons (Goo et al., 2006). Amini et al (1999), in their model of 
midbrain dopamine neurons, included CDI (steady state only) for N-type calcium 
but without any volt age- dependent inactivation. 

There have been a few voltage clamp studies of high threshold calcium currents 
in DRN SE neurons (Penington et al., 1991; Penington and Fox, 1995). Analysis of 
the I-V relations has been performed under the constant field assumption and the 
form of the current and its kinetic properties estimated. These analyses lead to the 
conclusion that the N-type current can be written 



Ca ; - Ca D e y o 1 

I N = k N m N h N V. (34) 

1 - e v o 

where m N and h N satisfy differential equations like (4) and (5). For activation 
variables we put 

mN >°° = i + e -(v-v Nl )/k Nl ( 35 ) 

and 

r mN = a N H N — r. (36) 

cosh((V - V N2 )/k N2 ) 

If the calcium N-type channels located on the soma (Penington et al., 1991) have 
similar properties to those examined in Joux et al. (2001), then they do not inacti- 
vate, as incorporated in the model of Komendantov et al. (2007). Some inactivation 
does occur however in DRN SE neurons and others (Cox and Dunlap, 1994; Lin et 
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al., 1997; Goo et al., 2006) so we put for inactivation variables 




1 



(37) 



I _|_ e (V-V N:i )/k N:i 



and a (large) voltage-independent time constant 



r hN = const. 



(38) 



Parameter values for I N 

Deducing kinetic parameters from voltage clamp data for In is also made difficult 
because of the many confounding factors such as temperature, mixtures of cur- 
rents, buffers and pipette solutions and nature of the charge carrier, which is rarely 
calcium. From the data of Figure 6 in Penington and Fox (1995), the following 
parameter values were estimated (voltages in mV, times in ms): = —13.5, 

k Nl = 9, a N = 0.305, b N = 2.29, V N2 = -20, k N2 = 20, V Ns = -50, k Nz = 20, and 
Th N = 1000. In their studies of variants of N-type channels expressed in sympathetic 
ganglia, Lin et al. (1997) found Vjvi values from -13.8 to -6.2, Vn 3 values from -53 to 
-43, r mN values from 3.3 to 5.8 and r hN values from 317 to 1567. In another study 
(unpublished data) we found V Nl = —15.2 and k Nl = 9. Furthermore, Tris-P0 4 
in the pipette shifted the activation by several mV to the left. These values may 
be compared with those employed by Komendantov et al. (2007), Vn 1 = — 11 and 
k^x = 4.2. From these estimates, the following set of parameters (Table 11) was 
chosen as standard for N-type. The value of k^ — 0.412 gives a peak current of 
about 3.8 nA with a clamp from a holding potential of -60 mV. 



Table 11: Standard set of parameters for I N 



Parameter Value Parameter Value 
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5 Calcium-dependent potassium currents 



Calcium entry into neurons (and other cell types) may activate certain calcium- 
dependent potassium ion channels which usually leads to hyperpolarizing effects. 
There are 4 main types of such channels in neurons, viz, the large conductance BK 
channel (or Kc a l.l) and the small conductance SK channels, SKI, SK2 and SK3 
(Kc a 2.1, Kc a 2.2 and Kc a 2.3, respectively). According to Camerino et al. (2007) 
BK channels are involved in a number of diseases including hypertension, coronary 
artery spasm, urinary incontinence, stroke, epilepsy and schizophrenia. 

Apart from the magnitudes of their conductances, the BK and SK types have 
differing activation properties and pharmacology. The activation of SK channels is 
purely calcium-dependent whereas that of BK channels depends both on calcium 
ion concentration and membrane potential. SK channels are selectively blocked by 
apamin from bee venom whereas BK channels are blocked by micromolar concen- 
trations of TEA and certain scorpion-derived toxins such as iberiotoxin. For reviews 
and discussions of the properties of these Ca 2+ -dependent potassium channels see 
Faber and Sah (2003), Muller et al. (2007), Fakler and Adelman (2008), Faber 
(2009) and Adelman et al. (2011). The role of BK channels is usually considered 
to be the hastening of repolarization after a spike with a consequent shortening of 
spike duration as demonstrated, for example, in rat CAl pyramidal cells (Shao et 
al.,1999; Gu et al., 2007; Loane et al., 2007), cerebellar Purkinje cells (Womack and 
Khodakhah, 2002), hippocampal granule cells (Muller et al., 2007, Jaffe et al., 2011) 
and dorsal root ganglion cells (Scholz et al., 1998). Blocking of BK channels may 
lead to bursting in some neurons (Traub et al., 2003). SK channels are involved 
in the medium duration (of order 100 ms) afterhyperpolarization following a spike 
which leads to a lengthening of the interspike interval. 

As noted above, calcium entry into neurons can be via several kinds of voltage- 
dependent calcium channel (VGCC) called L, N, P/Q, R and T-types, but not all 
are equally efficacious in the activation of Ca 2+ -dependent potassium currents. In 
CAl pyramidal cells Marrion and Tavalin (1998) showed that SK channels were 
only activated by L-type Ca 2+ currents, whereas BK channels were only activated 
by N-type currents, the latter also being the case for CAl granule cells (Muller et 
al., 2007). 

However, the situation differs from one neuron type to another. SK currents 
can also be coupled to N-type, P/Q-type, R-type or T-type (Sah and Davies, 2000; 
Faber and Sah, 2003; Loane et al., 2007). BK channels have also been found to 
be driven by R-type, L-type and P/Q-type (Fakler and Adelman, 2008). Sah and 
Davies (2000) contains a summary of couplings of various VGCCs with SK and 
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BK channels in various neurons. Close coupling of VGCCs to BK occurs when or 
because they are very close and there is evidence that the two kinds of channel 
are coassembled (Grunnet and Kaufmann, 2004; Fakler and Adelman, 2008). On 
the other hand, SK channels are not usually so tightly connected with VGCCs, 
allthough SK2 and T-type VGCCs are possibly co-assembled in some cells. In 
midbrain dopamine neurons, depolarization by L-type calcium currents occurs as 
part of the pacemaker cycle and the AHP in these cells is the result of the activation 
of SK channels by T-type calcium currents (Adelman et al., 2012). 

Calcium-dependent potassium currents have long been implicated as the basis 
of the long (several hundred ms) afterhyperpolarization (AHP) following spikes in 
DRN SE neurons and hence a major component of pacemaker- like activity in these 
cells (Vandermaelen and Aghajanian, 1983; Segal, 1985; Burlhis and Aghajanian, 
1987; Freedman and Aghajanian, 1987; Penington et al., 1992; Aghajanian and 
Sanders-Bush, 2002; Kirby et al., 2003; Beck et al., 2004). Direct experimental 
evidence has been provided by Freedman and Aghajanian (1987), Pan et al. (1994), 
Scuvee-Moreau et al. (2004), Rouchet et al. (2008), Crespi (2009) and Crawford et 
al. (2010). The recordings of Pan et al. (1994) and Scuvee-Moreau et al. (2004) 
show clearly the inhibiting effect of apamin on the afterhyperpolarization in DRN 
SE neurons. Since only SK channels are blocked by apamin, it is clear that the 
AHP is mainly due to SK channel activation. 

It is claimed that in rat brain BK channels are driven mainly by N-type VGCCs 
(Loane et al., 2007) so we will assume that this is the case in DRN SE neurons and 
that N-type calcium current contributes to the repolarization phase by activating 
BK channels. However, there is evidence for the colocalization and coassembly of 
L-type VGCC and BK in some parts of rat brain (Grunnet and Kaufmann, 2004), 
so despite the relatively small magnitude of the L-type Ca 2+ current in 5-HT cells 
of the dorsal raphe, there remains the possibility of an L-type-BK coupling. 

Anatomical investigations have been made of the distribution of BK and SK 
channels in mouse and rat brain, with similar findings in both species. In the first 
such study, Knaus et al. (1996) found that in rat, BK channel density was low in 
the brainstem relative to areas such as frontal cortex. A more detailed enumeration 
of BK densities in mouse brain (Sausbier et al, 2006) gives the numbers in the 
dorsal raphe nucleus as "few", compared with a maximum of "very high", thus 
corroborating the findings in rat brain. Stocker and Pedarzani (2000) found that 
in the rat DRN and MRN, SKI were absent, SK2 were at a "moderate" level and 
SK3 were "high". In B9 serotonergic neurons, both SK2 and SK3 were found to 
be "high". High levels of SK3 were found in the brainstem of mouse (Sailer et al., 
2004) and rat (Sailer et al., 2002). 
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Modeling the BK current 

The dual dependency of the BK current on Ca 2+ concentration and voltage makes its 
description a little more complex than that for purely voltage-dependent channels. 
Some of the first attempts at the inclusion of BK currents in spiking neuron models 
seem to be those of Yamada et al. (1989) for a cell of the bullfrog sympathetic 
ganglion, and Lytton and Sejnowski (1991) for cortical pyramidal cells. Schild et 
al. (1993) included a calcium-dependent potassium current to account for spike fre- 
quency adaptation in neurons of the rat nucleus tractus solitarii. More recently with 
additional data available, mathematical descriptions of calcium-dpendent potassium 
currents, including those through BK channels, have been handled in various ways 
(Traub et al., 2003; Xiao et al., 2004; Komendantov et al., 2007; Saarinen et al., 
2008; Jaffe et al., 2011). A simple model in which the calcium dependence is ignored 
was proposed by Tabak et al. (2011) and is explored in the simplified model in this 
article. 

The current through BK channels is designated by Ibk- The following equation 
describes this current on the assumption that the reversal potential is the Nernst 
potential for K + , 

Ibk = gBK,maxm BK {y - V K ), (39) 

where gsK.max is the maximal conductance and m BK is the activation variable. 
The constant field form may also be employed (Sun et al., 2004). It is assumed 
here that there is no explicit inactivation variable, although inactivation has been 
found when the BK channel has certain subunits (Xia et al, 1999) and has been 
documented in some other cases (Vergara et al. 1998; Sun et al., 2004; Gu et al., 
2007). The observation that the BK current often closely follows the calcium current 
indicates that the decline may be reasonably accurately described as deactivation. 
Furthermore, Marcantoni et al. (2010) found in mouse chromaffin cells that the L- 
type current through Ca v 1.3 channels is coupled to fast inactivating BK and as we 
have discussed above, it is most likely that in DRN SE neurons the L-type channels 
are .. 

With the simplification of no explicit inactivation 

drriBK _ m B K,oo - m B K 

dt DrriBK 

where the steady state activation and time constant depend on both V and Ca^. 

Empirically based estimates of the steady state activation function rriBK,oo(CcLi, V) 
have been available since Barrett et al. (1982). The data are usually characterized 
by fitting a family of Boltzmann functions of voltage at various values of Caj, where 
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(40) 



the half-activation potential and the slope factor depend on Caj. Many such data 
sets are available (for example, Cui et al., 1997; Scholz et al. 1998; Womack and 
Khodakhah, 2002; Sun et al., 2004; Thurm et al., 2005; Sweet and Cox, 2008) in 
a variety of preparations, which leads to very different parameter sets. Much of 
the data involves levels of Cai (and voltages) which are considerably higher than 
the ambient values expected in the normal functioning of DRN SE (or any other) 
neurons. This raises the question of how to interpret Ca, values in formulas for the 
activation functions and time constants. Since calcium ion concentrations just on 
the inside of BK channels, for example, just after the passage of Ca 2+ are much 
higher than the ambient averages (Fakler and Adelman, 2008), the latter need to 
be scaled up by a numerical factor cxbk when using formulas such as (41) and (44) 
below, which have been from data obtained using inside-out patch recording. 

On examination such data, those of Womack and Khodakhah (2002) were cho- 
sen, whereby a family of well-fitting Boltzmann functions was found to be, with 
half- activation potentials V B x,h and slope factors k B x, 

f V- v BK,h(Cai) \ - 1 

m BKt00 (C ai ,V) = il + e KsKicat) j ? ( 41 ) 

where 

V BK , h {C ai ) = -40 + 140e- (Ca *- a7)/7 (42) 

and 

k BK {Cai) = 11 + OMCtti, (43) 

with calcium ion concentrations in fiM and voltages in mV. The value of the half- 
activation potential V B K,h asymptotes to -40 mV at large calcium ion concentrations, 
displaying behavior similar to that shown in Latorre and Brauchi (2006). The steady 
state activation function described by these last three equations is plotted in Figure 
6. 

Empirical data on activation time constants for BK channels also exhibit much 
variability (Cui et al., 1997; Womack and Khodakhah, 2002; Sun et al., 2004; Thurm 
et al., 2005) which reflects the diversity of the preaparations. As a rare example 
where data at the same (Ca^, V) point are given in two studies, the time constant 
of activation reported at Ca, L = 10/zM and V = 20mV is 0.25 ms in rat inner hair 
cells (Thurm et al., 2005) and 1.3 ms in Xenopus motor nerve terminals (Sun et al., 
2004). Data are sometimes not given explicitly over ranges of values of Cai or V 
which would enable one to construct a function of these two variables suitable for 
use in a mathematical model. Since an action potential is only a few ms in duration 
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Figure 6: Steady state activation for BK channels as a function of membrane po- 
tential V and the Ca 2+ concentration in fiM. Adapted from data in Womack et al. 
(2002). 

it seems that for BK currents to affect this quantity the activation time constant 
should be not much greater than an ms at the corresponding values of Cdi and V . 
Consideration of the various data sets led to the use of that of Thurm et al. (2005) 
to approximate the activation time constant for BK, as the values were given over 
ranges of Cdi and V expected to occur during the action potential. Reasonable 
fits over the ranges of < Cai < 10 (in /xM) and —100 < V < 100 (in mV) were 
obtained on putting 

T mBK (Cai, V) = 1.345 - O.UCcii + V(0.0004Cai - 0.00455). (44) 

This yields a plane which is illustrated in Figure 7. 

Modeling the SK current 

This current is comparatively straightforward to describe, as its activation is pur- 
portedly only dependent on Cdi. We will assume that there is only one kind of 
SK channel, SK3, which was noted above to be prevalent in the DRN. If SK2 are 
present then they can be included in the same terms as their properties are very 
similar (Coetzee et al., 1999; Barfod et al., 2001). Thus 

ISK = 9SK,maxm S K(V ~ V K ), (45) 
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Figure 7: Activation time constant for BK channels as a function of membrane 
potential V and internal Ca 2+ concentration in juM. Adapted from Thurm et al. 
(2005). 



where 

dm SK m SKoo -m SK 

~ir = — r ' (46) 

though some authors have put m 2 SK instead of msK in Equ (45) (Yamada et al, 
1989; Destexhe et al., 1994). There is evidently no explicit inactivation process, 
as the decay of Isk is governed by Ca 2+ dynamics (Teagarden et al., 2008). The 
steady state activation is steeply dependent on Ccti and is written 

Ca n 

mSK >~ = CWtlQ (47) 

where n is the Hill coefficient and K c is the EC50 which is the value of Cai at which 
the activation msK,oo = 0.5. SK channels have K c values between 300 and 800 
nM and Hill coefficients between 2 and 5, where the ranges given by Vergara et al. 
(1998) and Stocker (2004) have been combined. The time constant T msK is about 5 
ms (Fakler and Adelman, 2008), the range 5-15 ms being given by Stocker (2004) 
and calcium dependency for this quantity was employed by Yamada et al. (1989) 
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and Destexhe et al. (1994). Sometimes the time constant has been ignored so that 
msK is always at its steady state value (Komendantov et al., 2007). Considering 
that time scales of several hundred milliseconds are relevant in the spontaneous 
activity of many DRN SE neurons, such an approximation would be reasonable 
except perhaps when the cells are in bursting mode. 



6 Hyperpolarization activated cation current, Ijj 

This current, which is is elicited by hyperpolarizations relative to rest, is slow to ac- 
tivate and does not inactivate (McCormick and Pape, 1990; McCormick and Hugue- 
nard, 1992; Pape, 1996; Robinson and Siegelbaum, 2003). In DRN SE neurons a 
similar current activating below -70 mV, which we denote by I H , was described by 
Williams et al. (1988). Ih may be enhanced by activation of certain 5-HT recep- 
tors thus preventing excessive hyperpolarization and tending to increase SE neuron 
firing rates (Aghajanian and Sanders-Bush, 2002). The current has been written in 
various forms and is often omitted in mathematical models, though it often plays a 
role in pacemaking - see the Discussion for details. 

In the absence of specific data we put, as in McCormick and Huguenard (1992), 

Ih = 9H,maxm H {V - V H ) (48) 

where V# is the reversal potential, usually taken as about -40 mV, toh satisfying 
an equation like (4). The steady state activation is given by 

1 

mH '°° ~ i + e (.v-v Hl )/k Hl > ( 49 ) 
with time constant well fitted by 

TmH = cosh((v- H v H2 ) /k H2 y (50) 



Parameter values for Ih 

The parameters for I H are slightly modified from those of Huguenard and Mc- 
Cormick (1992) to make the current activate at around -70 mV and have an ac- 
tivation time constant around 300 to 1000 ms (Williams et al., 1988). Thus the 
following set of parameters was chosen as standard for Ih- 
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Table 12: Standard set of parameters for J# 



Parameter 


Value 


Parameter 


Value 


v Hl 


-80 


V H2 


-80 




5 


k H2 


13 


a H 


900 


V H 


-45 



7 Fast transient sodium current, 

The only sodium current included is the transient Inu which, when blocked by TTX 
in DRN SE neurons, reduces spike amplitude by about 60 mV or more (Segal, 1985; 
Burlhis and Aghajanian, 1987). The current is given by the classical form 

iNa — 9Na,max m NahNa(y — V^a) (51) 

with activation variable rriN a and inactivation h^a- For the steady state activation 
we put 

mNa >°° = i + e -(v-v Nai) /k Nai ( 52 ) 
with corresponding time constant 

T m ,Na = a Na + b Na e~ ^ V ~ V ^ )/^- 2 ) ' ? (53) 

which fits well the forms used by some authors (McCormick and Huguenard, 1992; 
Traub et al., 2003) but not all. The steady state inactivation may be written 

hNa >°° = i + e (v~v Na3 )/k Na3 ( 54 ) 
with corresponding time constant fitted with 

r h , Na = c Na + d Na e-( (v - v ^ )/k ^y . (55) 



Parameter values for I^ a 

The properties chosen for fast sodium transient currents in neuronal models are 
diverse, especially with respect to time constants of activation and inactivation, 
which might be explained by the various types of Na v x channels, where x is 1.1 
to 1.9 (Catterall et al., 2005). Some examples which illustrate this are given Table 
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13. Note that the entry for the Ho dgkin- Huxley squid axon model gives potentials 
relative to zero as resting level, which was near -60 mV (Hodgkin and Katz, 1949; 
Tuckwell, 1988, Chapter 2). 



Table 13: Examples of parameters for I Na 



Source 










I'm, max 




Hodgkin (1952) 


+25 


10 


+2.7 


7.5 


0.5 


8 


Belluzzi (1991) 


-36 


7.2 


-53.2 


6.5 


0.3 


22 


Traub (2003) 


-34.5 


10 


-59.4 


10.7 


0.16 


0.16 


Komendantov (2007) 


-34.6 


6.2 


-61.6 


6.8 





29.1 


This study 


-33.1 


8 


-50.3 


6.5 


2 


8 



Voltage clamp experiments previously performed on a putative DRN SE neuron 
enable estimates of parameters for sodium current kinetics to be made. Current 
versus voltage values were obtained with a holding potential of -80 mV and steps 
to various test potentials from -60 mV to +60 mV. Using G/Gmax gave a half- 
activation potential (for m 3 ) of -30.73 mV and a slope factor of 4.62. However, it 
is useful to employ the methodology in Penington and Tuckwell (2012) to obtain 
estimates of the parameters of m NatOCJ and r m ^ Na as well as gNa,max- Current versus 
time curves for test potentials of -35 mV and -40 mV were obtained and that for 
-35 mV is shown in Figure 8. 

The maximum amplitude of the current and its time of occurrence are used to 
determine best fitting curves proportional to m 3 h as in Penington and Tuckwell 
(2012). This procedure yielded the following two equations (in nanoamps), 

1.12 = 62.73g Na , max m 3 Naj00 (-35) (56) 

0.657 = 75.33^ a , ma:E m 3 Va ,oo(-40). (57) 

It is convenient to assume that very nearly, mAr a ,oo(0) = 1, which yields, in con- 
junction with the experimental maximum current at the test potential V=0 and a 
supposed ratio of Th,Na to r m> Na at large depolarizations from two other experimental 
studies of sodium currents (Belluzzi et al., 1991; Yoshino et al., 1997), 

2.15 = 11.9^ O)mas m^ aiOO (0). (58) 

From this last equation the estimate is 5 , Ar a ,maa; = 0.18 /xS for a dissociated cell. It 
is then possible to estimate m7Va,oo( — 35) = 0.462 and m7v a ,oo(— 40) = 0.364. Using 
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Figure 8: Sodium current in a putative DRN SE neuron versus time for a test 
potential of -35 mV. 

the estimated values of m NajOQ at V=0, -35 and -40 mV, the Boltzmann parameters 
are obtained as in Table 14. 

To find r m jv a as a function of V we use the values estimated from voltage clamp 
data r m ^Na(— 35) = 0.15 ms and r m ^ a {— 40) = 0.2 ms. Considering the positions of 
the maximum of r m ^ a relative to the half-activation potential in other studies of 
sodium channel dynamics, it is reasonable to assume that the maximum of r m ^ a 
is in fact 2 ms at V = -40. Assuming too an asymptotic minimum of 0.05 ms at 
very large depolarizations and very large hyperpolarizations leads to the estimates 
of the four parameters in r m ^ a given in Table 14. 

For the steady state inactivation, in accordance with findings of Belluzzi et al. 
(1991), it is supposed that the half- inactivation potential is -50.3 mV at 17.2 mV in 
the hyperpolarizing direction from the half activation potential. The time constant 
for inactivation r h ^ Na was estimated from voltage clamp data at two test potentials 
to give 2.41 ms at -35 mV and 6.24 ms at -40 mV. Assuming as found in other 
studies that the maximum value of T^Na occurs at about 7 mV positive to the half- 
inactivation potential and that this maximum is 8 ms with a minimum of 0.5 ms 
yields the parameters for r h ^ Na and /iAr aj00 given in Table 14. Note that the value 
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of /cTVai has been reduced from its estimated value of 12.3 to 8, which is more in 
accordance with the average of the values found in other preparations. The standard 
set of parameters for I^ a kinetics is summarized in Table 14. 



Table 14: Standard set of parameters for I Na 



Parameter Value Parameter Value 



V Nai -33.1 V Na , -50.3 

k Nai 8 kjy a3 6.5 

a Na 0.05 c Na 0.5 

b Na 0.15 d Na 7.5 

V Na2 -40 V Na4 -43 

k Na2 7.85 k NaA 6.84 



8 Leak current, I Leak 

In the Ho dgkin- Huxley (1952) model, a leak current was inserted in the differential 
equation for V with its own equilibrium potential and conductance. This small 
current was stated to be composed of chloride and other ions, and the conductance 
did not depend on V. One motivation for including a leakage current is to take 
account of ion flows by active transport (pumps), although some authors take such 
transport into account explicitly. More recently, specific channels for potassium 
(Goldstein et al., 2001) and sodium (Tremblay et al., 2011) leak currents have been 
found. In the present model, as can be discerned from the properties of DRN SE 
neurons discussed above, there is not expected to usually be a zero ion flux at 
rest because many of these neurons fire spontaneously without afferent input. We 
therefore insert a token leak current which is carried by sodium and potassium ions. 
Based on a dynamical balance equation the following expressions are obtained for 
the leak current which is composed of a potassium ion and sodium ion contribution. 



Ileak Ix,leak ~\~ 1 Na,leak 



(59) 




(60) 
(61) 



(62) 
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9Na,leak — Tj 9K,leak (63) 

where Vk and Vno, are Nernst potentials, Vr is the resting potential and Ri n is the 
input resistance of the cell. 

9 Calcium ion dynamics 

The intracellular calcium ion concentration Ca,i varies in space and time throughout 
the cytoplasm. It may undergo increases due to the inward flow of Ca 2+ through 
VGCCs and because of release from intracellular stores, including endoplasmic retic- 
ulum and mitochondria, and release from buffers. Decreases occur by virtue of the 
pumping of Ca 2+ to the extracellular compartment, absorption by various buffers 
and the return of ions to the intracellular stores. Many of these processes are sub- 
jects of ongoing research so that accurate quantitative modeling is not able to be 
carried out with certainty. Calcium buffers have important effects on neuronal firing 
and synaptic transmission and altered expressions of them are associated with many 
diseases including Alzheimer's and Parkinson's, epilepsy, schizophrenia and depres- 
sion (Schwaller, 2007). Buffers can increase firing rate by reducing the amount 
of Ca 2+ available to activate BK and SK channels. We first discuss some of the 
relevant properties of buffers before formulating a quantitative description of Ca 2+ 
dynamics. 

The overall scheme for the rate of change in intracellular calcium ion concentra- 
tion has three components due to inward current Ic a through votage gated calcium 
channels, buffering and pumping, 

dCcii dCcii 
dt dt 

We will deal with each of these terms separately. 
Calcium buffering 

Calcium binding proteins are classified as either principally sensors, such as the ubiq- 
uitous calmodulin, or buffers, which bind Ca 2+ when sufficiently large local increases 
in Cdi occur. The main buffers in neurons are parvalbumin (PV), calbindin-D28k 
(CB) and calretinin (CR). Assuming that binding of buffer B to Ca 2+ is according 



+ 



dCa~ 



Bn f f 



dt 



+ 



dCd; 



dt 



(64) 



Pum/n 
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to the simple scheme 



Ca ++ + B — CaB, (65) 

k off 

with k on in units of /iM _1 s _1 and k ff in units of s _1 . Here and in the following, 
Ca ++ bound to B is denoted as simply Ca. The dissociation constant is 

K d =^ (66) 

in /iM and as pointed out by Schwaller (2010), its value usually far exceeds the 
resting level of Ca { so that most buffer is free in resting conditions. Fast buffers 
have a relatively large value of k on . Generally PV is considered to be relatively slow, 
CB is faster and CR is the fastest. Table 15 gives typical approximate (average) 
parameters for the three main buffers, based on several reports (Bortolozzi et al., 
2008; Cornelisse et al., 2007; Faas et al., 2007; Faas and Mody, 2011; Nagerl et al, 
2000; Schwaller, 2007, 2009, 2010). Here CBP is calcium binding protein. Note 
that various concentrations of other metallic ions such as Mg 2+ can in some cases 
alter the parameters for calcium binding. In the case of CR, the kinetics of binding 
depend on the state of occupancy and can change as Cai changes (Schwaller, 2009) 
which makes computational descriptions difficult. 



Table 15: Kinetic parameters for buffers 


CBP 


K d (nM) 


k 


(/iM-V 1 ) ^(s- 1 ) 


Remarks 


PV 


50 


19 


0.95 




CB 


195 


12 


2.34 


High affinity, slow 


CB 


490 


80 


39.2 


Low-medium affinity, fast 


CR 


68 


310 


21.1 


R sites 


CR 


28000 


1.8 


50.4 


T sites 


CR 


36000 


7.3 


262.8 


EF5 



The concentration of a buffer is an important variable and is unfortunately 
the least available, except for a few cell types such as hair cells (Bortolozzi et al., 
2008; Hackney et al., 2005) cerebellar neurons (Hackney et al., 2005 and Schwaller, 
2010, and references therein) and hippocampal neurons (e.g., Miiller et al., 2005). 
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Considering many reports and omitting extremely high values, the concentration 
ranges in neurons are given as approximately 50-150 //M for PV, 40-350 /iM for 
CB and 20-70 /iM for CR. Nonspecified buffer concentrations were given as 45 /iM 
(fast) and 250 /iM (slow) by Sah (1992). Computational models have often followed 
Yamada et al. (1989) who employed a buffer with a concentration of 30 /iM in the 
shell just interior to the membrane and 3 /iM elsewhere. The kinetic parameters 
were K d —1 /iM and k on =100 /iM _1 s _1 . Several other authors have followed the 
approach of Schild et al. (1993). 

Buffer types in the DRN 

With immunocytological methods, studies have been made of the gross details of the 
distribution of PV, CB and CR across various structures. These include rat brain 
(Celio, 1990; Resibois and Rogers, 1992; Rogers and Resibois, 1992), rat hindbrain 
(Arai et al., 1992) and macaque brainstem (Parvizi and Damasio, 2003). There have 
also been studies specifically targeting the DRN of the squirrel monkey (Charara 
and Parent, 1998) and chinchilla (Jaworska-Adamu and Szalak, 2009, 2010). In the 
results of most of these works, it is not possible to know with certainty which of the 
many kinds of cell in the DRN (see for example, Lowry et al., 2008) are the ones 
that contain the CBP. 

Certain generalizations have been made concerning the distributions of the three 
buffers under consideration. According to Resibois and Rogers (1992), PV is usually 
found in GABA-ergic neurons. Further, PV, CB and CR are mostly concentrated 
in different nuclei, but in some nuclei many neurons are positive for more than one 
of these three calcium buffer proteins. Also, CR is not generally found in major 
cell groups which are serotonergic such as occur in DRN. Schwaller (2007) claims 
that a greater variety of neurons express CB than PV and that, in accordance with 
Resibois and Rogers (1992) and Rogers and Resibois (1992), CB and CR are rarely 
found in the same cells. 

Concerning PV, (Celio, 1990) and Schwaller (2007) suggest little if any PV is 
found in the rat DRN and Parvizi and Damasio (2003) found no PV positive cells in 
the macaque mesencephalic raphe nuclei. However, inter-species differences prob- 
ably exist, because Charara and Parent (1998) found that in the squirrel monkey 
DRN, there were a few 5-HT-ergic neurons that contained PV. Jaworska-Adamu 
and Szalak (2009) found only weak PV immunostaining in the chinchilla DRN. It 
seems reasonable to conclude from the available evidence that PV is probably not 
frequently present in DRN SE neurons of the rat. 

In both squirrel monkey and macaque, many DRN neurons were found to contain 
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this CB. Jaworska-Adamu (2010) found intense staining for CB in nearly all DRN 
cells. It may be concluded that CB is a frequently occurring calcium buffer in 
the DRN, including the 5-HT containing cells. Concerning CR, it seems from the 
pictorial results of Charara and Parent (1998), that CR may be quite prevalent 
in the squirrel monkey DRN. In rat, Arai et al. (1991), found many small and 
medium-sized cells in the lateral parts of the DRN contained CR, with few in the 
central part. Jaworska-Adamu and Szalak (2009) also found CR containing cells in 
the chinchilla DRN, but it is not known if these were serotonergic. 

In summary, from the above observations, it may be concluded that in rat the 
DRN SE neurons are most likely to have CB as their main buffer, but that some 
cells may use CR. With modeling, one can test if the parameters for CB or CR lead 
to very different results if all other parameters are the same. 



Modeling of the buffering contributions 

The term ^p- appearing in Equ (64) will be described on the assumption that 

Buff 

the endogenous buffer is in fact calbindin-D28k. Many approaches to modeling a 
contribution such as this have appeared in recent years, some of which (Schmidt et 
al., 2003; Canepari and Vogt, 2006; Schmidt and Eilers, 2009) follow the relatively 
simple independent sites scheme set forth in Nagerl et al. (2000) and some, such 
as Nadkarni et al. (2010) and Modchang et al. (2010) use more complete and 
hence complicated kinetic schemes for CB to Ca 2+ binding. Pursuing the simpler 
description, we let the concentration of the k-th buffer site, with k = 1,2, 3, 4, be 
Bkit) at time t and let the reaction with Ca 2+ be 



Ca ++ + B k ^CaB k , 



(67) 



where now the more succinct symbols f k and b k are used for forward (on) and 
backward (off) reaction rates. Applying standard chemical kinetics formulas gives 
the contribution to the rate of change of intracellular calcium ion concentration 
C ai (t) 

(,( '"' % * b k [CaB k ](t) - f k B k (t)Coi(t) 



dt 



l Buff 



k=l 



(68) 



where [CaB k ](t) is the concentration of bound k-th buffer. Following Yamada et 
al. (1989) we let the total (bound and free) concentration of buffer site k be the 
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(assumed) constant B ktot and note that 

B k , tot = B k {t) + [CaB k ](t), 



(69) 



which enables one to eliminate the bound buffer-calcium terms. The buffering 
contribution becomes 



dCdi 
dt 



= E 

Buff k=l 



b k (B k;tat - B k (t)) - f k B k (t)Coi(t) 



The 4 subsidiary equations are 
dB k 



dt 



b k [B kM -B k (t)] - f k B k (t)Ca t (t). 



(70) 



(71) 



These previous two equations enable one to march forward in time, for example 
with an Euler scheme, in which, the buffering contribution only gives 

4 



Ca t ((j + l)Af) = Ca t (jAt) + At ^ 
and 



b k [B kytot - B k (jAt)] - f k B k (jAt)Ca t ( 3 At) 

(72) 

B k ((j + 1) At) = b k [B Ktot - B k (jAt)] - f k B k ( 3 At)Ca, t ((jAt). (73) 



Calcium currents 



The term ^ 

at 



ICa 



in (64) has several components. The presence of L-type, N-type 



and T-type Ca + currents has been verified in DRN SE neurons. A further calcium 
current, called Ic a ,R ap he was also found with similar electrophysiology to the N-type 
current, but it was not blocked by u- Conotoxin-GVIA (Penington and Fox, 1995). 
In the quantitative treatment of calcium currents we treat N-type and Ic a ,Ra P he all 
as N-type. 

Because the calcium-activated potassium currents may distinguish the type of 
Ca 2+ channel which activates them, we keep track of the changes in the various 
calcium currents which have a total designated by 



lea — II + In + It, 



(74) 



and must also decompose the intracellular calcium concentration into its sources. 
It is assumed that the resting level Ca« r is the value at t — and that 



C ai (t) = Ca i;R + Ca^ L (t) + Ca hN (t) + Ca hT {t). 



(75) 



45 



However, as a simplifying assumption the total intracellular calcium concentration 
is used if the constant field expressions are employed and in determining the rates 
of buffering and pumping. 

The usual and often-used assumption is made that the change of the calcium 
concentration due to flux through the voltage-gated channels occurs within a thin 
shell just inside the cell membrane with ascribed volume v (Standen and Stanfield, 
1982; Yamada et al., 1989), which might depend on the type of current but this 
complication will not be included. Then,with F denoting Faraday's constant, for 
the rate of change in internal calcium due to the L-type current 



dCcii^L 



dt 2Fv 



(76) 



and similarly for the N-type and T-type currents. A calcium source factor CSF with 
a maximum value of 1 and a minimum value of is multiplied by the current terms 
to limit the magnitude of the internal calcium ion concentration, if necessary. 



Calcium pumping 

Intracellular Ca 2+ levels are kept low after surges of current through VGCCs and the 
maintenance of such levels is brought about by buffering and active transport to in- 
tracellular stores or to the extracellular space. The two main pumping mechanisms 
are the PMCA (plasma membrane Ca 2+ ATP-ase) pump and the NaX Na + -Ca 2+ 
exchange pump (Carafoli, 2002). Only a few models of single neurons (for exam- 
ple, Lytton and Sefnowski, 1991; Amini et al., 1999) have included the NaX pump 
but all include versions of the PMCA pump. Some authors have employed a lin- 
ear form (Komendantov et al, 2007; Cornelisse et al., 2007) or voltage-dependent 
forms (Rybak et al., 1997) for the latter but most often the form utilized is the 
simple Michaelis-Menten one (for example Good and Murphy, 1996; Lumpkin and 
Hudspeth, 1998) 



dCdi 
dt 



Pump 



where K s determines the maximal pumping rate and K m is the dissociation constant 
which is the value of Caj for half-maximal rate. It seems that the reasoning of 
Lumpkin and Hudspeth (1998) for using this simple form and not including NaX, 
is sound, namely that there is a lack of reliable data for the required parameters in 
more complex models. 
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10 Simplified model 



It is convenient to first see if some of the firing properties of DRN SE neurons can 
be predicted with a model which has the main elements of the model described 
in Sections 2 to 9 but has a few simplifications. In the first instance, the latter 
consist of (a) ignoring the BK current (b) for calcium currents, using the form of 
(18) with an assumed fixed calcium reversal potential rather than the constant field 
form (c) omitting the M-type potassium current (d) simplifying the equations used 
to describe the buffering of calcium to those employed by Rybak et al. (1997), so 
that the differential equation for internal calcium concentration becomes 

*** ^sFv l + i s )}-j B ®-K. Ca ' 



dt v ' 2Fv Cai + Kr, 

where the fraction of calcium which is bound is 

B t ot 



PB(t) = 



Ca,i + B tot + K d ' 



Kd being the dissociation constant defined in (66). In addition, an applied current 
/i, which may be positive (hyperpolarizing) or negative (depolarizing) is added. This 
applied curent can be experimentally applied or have its origin in synaptic inputs 
such as the noradgrenergic input from locus coeruleus or inhibitory inputs from 
neigbouring DRN SE cells or within nucleus GABA-ergic neurons. 



Equation (6) for the voltage now becomes 
dV 

C— = -[I A + I KDR + I T + I L + I N + I H + I Na + I SK + I Leak + lA (78) 

where the equations describing the individual components are given in the individual 
sections describing them. For time constants of activation and inactivation variables, 
a further simplification is sometimes made that they are not voltage dependent 
as given in their text descriptions. For example, the time constant r mi jv Q , which 
depends on 4 parameters in Equ. (53), may be replaced with the single parameter 

(79) 

etc. 
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Table 16: Activation and inactivation parameters 
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Table 17: Cell properties, Ca dynamics, Equilibrium potentials 



V K 
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Ca Q 
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c 
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A 
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0.7 


d 
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n 






In the simplified model there are 89 parameters, many of which have been given 
in the text. Tables 16-18 give the set of parameter values, referred to as Set 1, which 
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Table 18: Maximal conductances 



9Na,max 2 9N,max 0.5 

9KDR,max 0.5 Q A,max 0.479 

9T,max 0.0825 9lH,max 0.005 

9L,max 0-0825 g S K,max 0-01 



are first employed in the simplified model. This is a trial set which will be modified 
as necessary in the sequel. Table 16 contains the parameters required to describe 
the steady state activation and inactivation functions as well as the time constants 
for the various channels. Table 17 gives ionic equilibrium potentials and several cell 
properties including parameters required to describe calcium buffering and pump- 
ing. Finally Table 18 gives 8 maximal conductances, being for the 7 voltage-gated 
channels and the calcium-activated potassium current (SK). All potentials are in 
mV, all times are in ms, all ionic concentrations are in mM, all conductances are in 
fiS, A is in square microns, d is in microns, C is in nF, and F is in .coulombs/mole. 

40 
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Figure 9: Isolated spike in the simplified model with set 1 parameters as in Tables 16, 
17 and 18, and a depolarizing applied current fi = —0.17 nA. The resting potential 
is -60 mV. 



49 



11 Computed solutions 



The differential equations of the model were integrated with an Euler method of 
step size 0.004 ms, the results being indistinguishable from those obtained by Runge- 
Kutte techniques. Initial conditions were chosen as values at resting potential and 
resting ion concentrations. In the following the equilibrium potential for Ca 2+ was 
set at +60mV. 

Preliminary results for the 9-component model (78) 

With no applied current the results for Set 1 parameters showed no evidence of 
spontaneous spiking. The membrane potential drifted down from resting value of 
-60 mV to almost be at about -65 mV in 25 ms. It is clear that the depolarizing 
forces of low threshold calcium It and fast sodium current were overwhelmed by 
the negative drive from the transient A-type potassium current I a and the delayed 
rectifier potassium current. With an applied depolarizing current [i < 0, the same 
situation prevailed until it reached a level of ji = —0.17 nA. With this value of //, 
as depicted in Figure 9, a spike was emitted at about t=55 ms, with maximum 
value +27.5 mV, minimum -80 mV and a duration D (always taken to be the spike 
width at -40 mV) of 1.4 ms. The magnitudes of both the fast sodium current 
and the delayed rectifier potassium current were both very large (about 18 nA). 
Even after a very long time interval, there was no subsequent spike as V tended 
to remain at a new equilibrium value around -71 mV. There was a pronounced 
afterhyperpolarization as found in DRN SE cells, but the behavior of V after the 
spike was not typical for these neurons. 

Consideration of simplest model with only Na and K 

With a multi-component neuronal model there are of course a large number of vary- 
ing properties for different choices of parameters. Just using three values of each 
of 89 parameters would result in over 700,000 combinations. With so many param- 
eters it is difficult to see which are responsible for important spiking properties. 
Many chosen combinations led to satisfactory spike trains with properties similar to 
experiment, but often with durations which were unacceptably long. Other combi- 
nations led to spontaneous activity but spikes had uncharacteristic large notches on 
the repolarization phase. Doublets and triplets were often observed, and although 
these are found in experimental spike trains, it was desirable to find solutions depict- 
ing the regular spiking and relatively smooth voltage trajectories usually observed. 
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Hence a gradual increase in complexity was employed starting with the usual basic 
two currents in action potential generation. 

In this subsection attention is focused on the spike properties of duration and 
to a lesser extent inter-spike interval (ISI). In subsequent sections the various other 
component currents are added. We illustrate with two choices of parameters for the 
sodium and potassium variables m Nat00 , hjy ajOQ and n^, one set being chosen from 
the above parameter set and another based in part on those of Belluzi and Sacchi 
(1991) which we denote by Set 2. For the latter we also use the resting potential 
from Kirby et al. (2003) and the cell capacitance, both being for rat DRN SE cells. 
The two sets of parameters are summarized in Table 19. 

Table 19: Two basic parameter sets for the Na-K system 



Parameter Set 1 Set 2 
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Both of these reduced model parameter sets led to repetitive spiking with the 
addition of a small depolarizing current. Typical spike trains are shown in Figure 
10 and Table 20 lists some of the details of the spike and spike train properties. 
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Figure 10: Repetitive spiking in the reduced model (Na-K) for the parameter sets 
of Table 19, Set 1 (top part) and Set 2 (lower part). In both cases /x = —0.05. 

Spiking for the second parameter set has a lower threshold for (repetitive) spiking, 
a longer ISI at threshold, a longer spike duration and a larger spike amplitude. For 
both parameter sets, most of these spike properties are in the ranges observed for 
DRN SE neurons so that such sets could form the basis of a more complete model 
in accordance with Equ. (6) as explored below. 

Figure 10 shows well defined spikes with abruptly falling repolarization phases 
to a pronounced level of hyperpolarization followed by a steady increase in depolar- 
ization until an apparent spike threshold is reached. DRN SE neurons are usually 
characterized as having long plateau-like phases in the latter part of the ISI and 
this is a feature shown in Figure 11 for spikes elicited near the threshold for spiking 
for both sets of parameters. The plateau for the second set is nearly three times as 
long as that for the first set. 

Graphs of the frequency of repetitive spiking versus depolarizing input current 
are shown in Figure 12. In both cases it seems that at a particular value of /x the 
frequency jumps from zero to values of 3.0 Hz for set 1 and 1.1 Hz for set 2, rather 
than increasing continuously from zero. Thus these models with the chosen param- 
eters woukd be classified as type 2 neurons (Hodgkin, 1948; Tateno et al., 2004) 
whereby the instability of the rest point is due to an Andronov-Hopf bifurcation 
rather than a saddle-node bifurcation as in Type 1 neurons. These results can be 



52 





Set1:n=-0.0342 






Plateau -260 ms 





200 250 300 350 400 450 500 550 600 650 



V(mV) 



20 













Set 2: |i=-0.018 




-20 








-40 




-A 




-60 












Plateau -700 ms 




-80 









700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 



Figure 11: Showing long plateaux in repetitive spiking in the reduced model (Na-K) 
for the parameter sets of Table 19, Set 1 (top part) and Set 2 (lower part). In both 
cases n is very close to the threshold for firing. 



compared with those for the complete simplified model in Section 12. 



Table 20: Spike train properties 



Property Set 1 Basic set 2 

Spike threshold \i = -0.0342 Ji = -0.018 

ISI at threshold 331 ms 948 ms 

Spike duration (-40 mV) 1.6 ms 2.9 ms 

Max V +8 +19.4 

Min V -90.0 -91.2 



Four- component model with a high threshold Ca 2+ current 

and I S K 

We wish here to include in the model the spike-induced influx of calcium ions and 
calcium-activated potassium currents. For the former, we reduce the components 
described in Section 4 to just N-type Ca 2+ and for the latter only an SK current 
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Figure 12: Frequency of action potentials versus magnitude of depolarizing current 
for the two above parameter sets with only two components. Note the different 
ordinate scales. 

is included along with calcium buffering and pumping. Thus these results would 
apply when the currents I a, It, Ih are blocked. For these computations three of 
the parameters of the N-type Ca 2+ activation and inactivation were changed to 
Vat 1 = —20, Vn 3 = —50 and /cjv 3 = 10, representing a general high threshold Ca 2+ 
current. 

Two sets of results for V, the internal Ca 2+ concentration, and Isk are shown in 
Figure 13, where also the near threshold value of \x = —0.035 was employed. In the 
results shown in the left column of Figure 13, the following maximal conductances 
were employed g N ,max=0-±, gNa,max=2, gKDR,max=0-5, gsK,max=0-003, so that the 
sodium and potassium conductances are the same as in Set 1 of Table 19. The ISI 
is about 550 ms, the maximal values of Ijv a , Ikdr, Ikca and In are about 10, 13, 
0.03 and 0.7 nA, respectively and the maximal intracellular Ca 2+ concentration is 
only about 400 nM, the resting level being 50 nM. There is a plateau in the latter 
phase of the ISI as the threshold for the subsequent spike is approached gradually. 

In order to reduce the magnitudes of ijva and Ikdr, the following set of maximal 
conductances was employed: gN,max=0-l (as before) gNa,max=l, gKDR,max=0.25, 
5 , sii",maa;=0.015,and the results are shown in the right-hand column of Figure 13. 
Whereas the ISI is little changed with a value of about 544 ms, the membrane 
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potential during the ISI has a decidedly different form in that it remains at the most 
hyperpolarized levels for more than half of the ISI and then turns fairly sharply 
upward to reach threshold for firing at similar values (-50 mV) as for the first 
set of conductances. Now however, the maximal values of In a, Ikdr, Isk and 
In are about 3.5, 5.5, 0.1 and 0.9 nA, respectively and the maximal intracellular 
Ca 2+ concentration is about 600 nM. The form of Isk is also quite different in the 
two cases considered, there being a rapid increase (in absolute value) and equally 
rapid decrease around the time of occurrence of the spike for the second set of 
conductances. 
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Figure 13: Showing the following three quantities versus time during spiking : volt- 
age (V) (top), intracellular calcium ion concentration (Ca^) (middle) and calcium- 
activated potassium current (Irca) (bottom). In the left-hand column the key 
parameters are g N ,max=0.l, g N a,max=2, gKDR,max=0.5, gsK,max=0.003 (set 1 above) 
and in the right-hand column, the same conductance for N-type Ca 2+ but gNa,max=^, 

gKDR,raax=0.25, gSK,max=0M5 (set 2). 



Changes of calcium pump strength 

If the calcium pump is reduced in strength so that the intracellular calcium con- 
centration does not decline so rapidly after Ca 2+ influx from the high threshold 
current In, then the value of msK,oo in (46) remains high and hence Irca is slow 
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to decline. This leads to a prolongation of the ISI as is seen in the records shown 
in Figures 14 and 15. In Figure 14, which shows membrane potential versus time 
during the ISI, the top two records were obtained with the first set of conductances 
and the bottom two records with the second set, as described in the previous two 
paragraphs. It can be seen that despite the lengthening or shortening of the ISI by 
decreasing or increasing the calcium pump strength, respectively, the general form 
of the response is not altered. That is to say, for the first conductance set there 
is still a prolonged and gradual plateau-like approach to threshold, whereas for the 
second set, V remains near its most hyperpolarized levels for half or more of the ISI 
and swings upward to threshold rather sharply towards the end. The corresponding 
time courses of the intracellular Ca 2+ concentration are shown in Figure 15. 
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Figure 14: Effects of increased and decreased calcium pump strength, described by 
the parameter K s in Equ. (77) on membrane potential during spiking. In the top 
two records gNa,max, gKDR,max, 9sK,max , In conductances from set 1 are used and in 
the bottom two records the set 2 values as described in the caption of the previous 
figure. The numbers on the plots indicate increases (1.5) or decreases (0.5) in K s 
values relative to the standard value. 

Smaller I^ a and Ikdr and the AHP 

In the above simulations for set 1 and set 2 parameters, the maximal sodium con- 
ductances were 2.0 /xS and 1.0 /jS respectively and the maximal values of gKDR,max 
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Figure 15: Intracellular Ca 2+ concentrations corresponding to the voltage trajecto- 
ries in the previous figure for two sets of conductances and two pump strengths. 

were 0.5 and 0.25 /iS, respectively. Although the properties of the spikes as shown 
in Figures 10, 11 and 14, 15 could serve as a basis for observed spikes, the sodium 
and potassium currents were fairy large, being of the order of 10 nA. 

A sodium current of about 9 nA was given for hypothalamic magnocellular neu- 
roendocrine cells (Komendantov et al. 2007), but judging by measured or computed 
currents in cells expected to be similar to DRN serotonergic cells (Milescu et al., 
2008; Milescu et al., 2010), the sodium currents are likely to be of order 2-3 nA. 
Further suppport for this estimate comes from analysis of a clear image of a spike 
in Kirby et al. (2003). The derivative of the voltage trajectory at the upswing of 
the spike is estimated at about 63 volts per second. From the same report the given 
resistance and time constant lead to an estimate of 0.0886 nF for the capacitance, 
which combined with the value of dV/dt give a current of about 5.6 nA. Assuming 
this is about 90% sodium current results in an estimate of 5 nA. On the other hand, 
since the capacitance estimate may be high (Tuckwell, 2012b), if the value 0.04 
nF is used then the estimate of the sodium current is about 2.4 nA, in accordance 
with the estimate from the articles of Milescu et al. (2008, 2010). In conclusion, 
the maximal sodium current in DRN serotonergic neurons is likely to be about 2-3 
nA with the possibility of values a high as 5 nA. One expects the delayed rectifier 
potassium current to be of approximately the same or smaller magnitude. 
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In fact, recall the value gNa,max=^-^8 /xS estimated in section 7 which was for a 
dissociated cell at room temperature. The dissociated cell has a small capacitance 
of perhaps only 20 pF (Penington and Fox, 1995) which means a capacitance of 9 nS 
per pF at room temperature or perhaps about 13.5 nS per pF at body temperature. 
Thus an average sized cell of capacitance 40 pF (set 1 cell) would have (?va,max=0.54 
/xS and a larger cell with a capacitance of 88.61 pF (set 2 cell, Table 19) a value 
9Na,max = ^--^ A*S. The first of these values is considerably less than the ad hoc value 
of 2 /xS in Tables 18 and 19, but the second value is close to the value in set 2 of 
Table 19. In a similar vein, the conductance for delayed rectifier potassium current 
was given as 0.641 nS/pF in Section 3.3 which translates to 0.0256 /xS for an average 
cell with C = 40 pF and 0.0567 /xS for a larger cell with C=88.61 pF. 

In investigating the effects of changing various parameters on the spike train 
properties, a systematic approach would require a large amount of space, as each 
of the about 90 parameters does have an important influence. Small changes in 
just one parameter can radically alter the nature of spiking. With the above new 
estimates for gNa,max and gKDR,max, a number of runs, reported in Table 21, were 
performed. 

Runs 1-5 of Table 21 

With the estimated values for QNa,max (0.54) and QKDR,max (0.0256), model responses 
(with just Na and K) to a depolarizing current of 0.05 nA were computed using the 
remaining Set 1 values of Table 19 for sodium and potassium DR currents. This gave 
rise to a rapidly oscillating solution, similar to action potentials at the beginning, 
but dying out after 20 ms at an equilibrium value of about -20 mV. Changing the 
value of nfc to 2, which moves the Irdr activation in the depolarizing direction, led 
to an even briefer oscillatory response, stabilizing again at about -20 mV. With nk 
again at 1, the time constant for potassium activation (17) was increased by putting 
&KDR = 14 so it had a maximum of 15 ms. Then with /x still at -.05 nA, repetitive 
spiking occurred with an ISI of 15 ms, duration (-40) of 5 ms, V max = 11.5 mV, 
Vmin = -55 mV, I Na>ma x = 2.0 nA and I K DR,max = 1-4 nA. 

Several sets of results are summarized in Table 21. In all cases, n\. — 1, /x = 
—0.03 nA and the value of bxDR = 14 in the potassium activation time constant. 
In all runs except the 8th, a regular train of action potentials emerged. In Run 
1, as in the case with /x = —.05 nA, the time constants for sodium activation and 
inactivation (called just r m and r h here) are constant (Table 19, 1st column). The 
ISI is slightly longer at 19 ms, and generally there are only small changes in spike 
train properties, including the maximum sodium and potassium currents (called just 
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INA and IK in Table 21). The duration D is the same at 5.0 ms. In Run 2, the time 
constants for sodium are given their full voltage-dependent forms (denoted by vd in 
Table 19) as in Equs (53) and (55) with parameters as in Table 14. This increases 
Ino, significantly, leads to more extreme values of V, shortens the duration to D=3.9 
ms and increases the ISI to 22 ms. With the maximal sodium conductance gNa,max 
reduced by 25% to 0.41 fiS, (Run 3), the ISI is hardly changed, but the maximal 
voltage and sodium current are significantly less. The duration is reduced by 0.1 
ms to 3.8 ms. A further reduction in gNa,max to 0.35 /xS (Run 4) leads to a sodium 
current about as estimated from experiment at 5.4 nA, without much change in the 
remaining properties in Table 19. Run 5 also includes the leak current (see section 
8), all other parameters being as in Run 4. The effect on the ISI is large as it 
increases from 23 to 103 ms, with a concomitant drop in maximal sodium current. 
The remaining spike properties are not greatly different, except D is greater at 4.6 
ms. 



Table 21: Effects of changing some parameters on spiking 



Run 


GNA 


GK 




Th 


Leak 


GN 


GSK 


Vrnax 




ISI 


D 


INA 


IK 




1 


0.54 


0.0256 


C 


c 








10 


-55 


19 


5.0 


1.8 


1.4 




2 




» 


vd 


vd 








30 


-60 


22 


3.9 


7.4 


1.9 




3 


0.41 














24 


-60 


22 


3.8 


5.4 


1.7 




4 


0.35 


ii 












20 


-58 


23 


4.0 


4.6 


1.5 




5 


n 


ii 






yes 






7.4 


-56 


103 


4.6 


2.9 


1.0 




6 


n 










0.003 


0.012 


7.4 


-56 


255 


4.8 


2.9 


1.0 


100 


7 


n 


ii 








0.006 




7.5 


-65 


491 


4.7 


2.9 


1.05 


170 


8 


5) 


0.032 








0.006 






No 


AP 










9 


0.41 


ii 








ii 




11 


-64 


383 


4.1 


3.6 


1.4 


156 


10 


n 


ii 


C 


c 




ii 




-13 


-79 


750 


10 


0.65 


0.33 


300 


11 


0.35 


0.0256 


vd 


vd 




0.012 




-7.5 


-81.5 


1100 


5.3 


2.8 


1.05 


393 


12 


1.62 


0.0768 


vd 


vd 




ii 




39 


-79 


475 


2.4 


17 


4.8 


295 


13 


1.08 


0.0768 


vd 


vd 




ii 




33 


-76 


422 


2.4 


12 


4.3 


248 



Runs 6-13 of Table 21 

In Run 6, the high threshold calcium current 1^ is present with a maximum (de- 
noted by GN in the Table) of gAr im a:r=0.003 /iS, and at the same time the calcium- 
dependent potassium (SK) current is switched on with a maximal conductance 
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(denoted in the table by GSK) of gsK,max= 0.012 //S. Without any change in the 
resulting sodium and potassium delayed rectifier current, nor in the maximum and 
minimum votages, the ISI is more than doubled to 255 ms, the spike dration D is 
slightly increased to 4.8 ms and the internal calcium concentration rises to 100 nM 
at its maxima. 

In Run 7, the only change in parameters from Run 6 is the doubling of the N-type 
calcium conductance to gN,max= 0.006 fj,S, with a consequent deeper AHP, a slightly 
less duration D, and a much longer ISI of 491 ms. The maximum internal calcium 
concentration is now 170 nM. The sodium and delayed rectifier potassium currents 
are little altered. The spikes have an appearance similar to those of some DRN SE 
cells, and some of their properties are shown in Figure 16. The properties shown are 
V, the internal calcium concentration, the fast sodium and N-type calcium currents 
and the potassium currents Irdr and Isk- 

Because the AHP was not as deep as usually found in DRN SE cell spikes, the 
value of gKDR,max was increased by 25 % in Run 8, with all other parameters the 
same. The result was that no spiking at all emerged. To counter this, in Run 9, 
the value of gNa,max was increased back up to 0.41 as in Run 3. Spiking was re- 
established but compared to the Run 7 spike, the maximum of V was increased to 7 
mV and the minimum increased to -64 mV, whereas the spike duration was reduced 
to 4.1 ms and the ISI was also decreased to 383 ms, the maximum Cai was only 156 
nM, and the maxima of I^ a and Irdr were greater at 3.6 and 1.4 nA, respectively. 

Another run (Run 10) with the same parameters was performed except with 
constant rather than V-dependent time constants for I^ a , which resulted in a much 
greater ISI at 750 ms and a much greater spike duration of 10 ms. At the same 
time Ino, and Irdr were very small. Whereas the minimum voltage was -79 mV, 
the maximum was also decreased to -11 mV. Run 11 is the same as Run 7 but 
with increased gN tmax — 0.012. This causes a lengthening of the ISI to 1100 ms, the 
duration is somewhat larger at 5.3 ms and the internal Ca 2+ increases to 393 nM, 
being more than doubled. 

Effects of changes to parameters for I N 

In order to determine the influence of the parameters for activation and inactivation 
of the N-type calcium current, the following (denoted by N-type set 2) were used 
(with all else as in Run 11), as estimated from a voltage-clamp study (Penington 
and Fox, 1995): V Nl = -13.5, k Nl = 9, a N = 0.305 , b N = 2.29, V Ni = -20, 
fcjv 2 — 20, Vn 3 = —50, /ctv 3 = 20 r^jv = 1000. The data are for barium as charge 
carrier, room temperature and with EGTA and may be different for native cells at 
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body temperature. 

For the results, defining the summarizing 7-vector 

X W^maxi ^mint ^ i ISI, I N a,maxi lKDR,maxi C Climax] , 

this gave, X = [87.3, —65.6, 5.1, 500, 2.85, —1.05, 175], differing significantly only in 
Vmin, ISI and Cai^ max . In this and all of the above runs 1-11, the spike duration was 
longer than that expected for these cells, in contrast with the value 1.6 ms in Table 

20, Set 1. 

With the time constants and (larger) conductances for I Na and Irdr as in set 
1 of Table 19, and with N-type set 2, a more satisfactory result was obtained, 
X = [26.8,-80,3.1,587,6.3,-5.1,349]. The results of these last two runs are not 
given in Table 21. Similarly with the original set 1 for N-type, Run 12 in Table 

21, with increased gNa,max and gKDR,max, a very reasonable duration D=2.4 ms was 
obtained but with very large I^ a and a smaller ISI. With gNa^max reduced by 33%, 
the duration was unchanged - see Run 13 in Table 21. With \x = in the last run, 
no spikes emerged, nor when the value of gKDR,max was reduced to 0.0512 or further 
to 0.0384. 
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Figure 16: Spike properties for Run 7 (see Table 21) described in the text and with 
parameters as given in Table 19. The ISI is 491 ms, or a frequency of 2.04 Hz. 
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Effects of Ibk and changing capacitance and area 

In this subsection the effects of altering the cell capacitance C and area A are 
briefly investigated as well as the effects of including a simplified model for BK 
current (Tabak et al., 2011). The results are given in Table 22, being labeled 
SI to S6. In all these runs the following maximal conductances are employed: 
9Na,max=0.675, g K DR,m a x=0 .0768, 0jv, mO x=O.O24 (double that in the last line of Table 
21) and g SK , m ax=0M2. 

In SI, the dramatic effect of increased gN,max is seen, especially in that the ISI 
is much longer at 1253 ms, well within the range of firing rates of many DRN SE 
neurons. In S2, the cell capacitance is increased by 25% and this change induces an 
irregular pattern of ISIs, at least for the first several seconds, the pattern being long, 
short, long, short, with values given in column 7 of Table 22. The spike duration is 
also increased as is also the peak internal Ca 2+ concentration. There appeared some 
small undulations in the potential near the end of the ISIs as spikes were approached. 
When the total area A is increased by 25%, as in S3, the train was regular again, 
with a decreased ISI, a little-changed spike duration and a slightly smaller internal 
Ca 2+ concentration. In a case of increases in both area and capacitance, (S4), a 
regular train arose with an ISI of 932 ms and spike duration 3 ms. 

As the spike durations in most of the runs so far have been at the long end 
of the range of those commonly reported for DRN SE neurons, it is of interest to 
consider the effects of potassium BK currents, which are often posited to play a role 
in speeding up the process of spike repolarization. In the mathematical model for 
BK currents presented in Section 5.1, the values of internal Ca 2+ required to make 
the conductance significant are above those obtained in the computations reported 
above. Hence a simpler more direct model was employed as in Tabak et al. (2011) 
whereby the steady state activation depends only on voltage. The half-activation 
voltage was set at -20 mV, the slope constant at 2 mV, both as in Tabak et al. 
(2011) and the time constant at 2 ms, regardless of voltage, which is at the low end 
of the range in Tabak et al. (2011). There was, as usual, no inactivation for BK 
channels (see Section 5.1). 

Runs S5 and S6 include BK current with conductances equal to 1/3 and 2/3, 
respectively, of the potassium delayed rectifier current. The effects of the inclusion 
of BK current here were surprising as even though the spike duration was reduced 
from 2.65 ms to 2.4 and 2.2 ms for S5 and S6 respectively, the spike train became 
irregular. For the smaller BK conductance (S5) , long and short ISIs alternated, 
being 1147 and 876 ms respectively and the peak internal Ca 2+ concentration also 
alternated. For the larger BK conductance (S6), the pattern of ISIs was long, short, 
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short, long, short, short etc. at least for the first several seconds. 



Table 22: Induced firing with \x = —0.03: changes in C,A and effects of BK 



Run 


GBK 


C 


A 


V m ax 


V m in 


ISI 


D 


INA 


IK 


Ctti^max 


SI 




0.04 


4000 


20.7 


-81 


1253 


2.65 


6.9 


-3.6 


380 


S2 




0.05 


4000 


16.7 


-81 


1360, 1678 


3.0 


6.6 


-3.4 


400 


S3 




0.04 


5000 


21.2 


-80.5 


826 


2.7 


7.0 


-3.7 


344 


S4 




0.05 


5000 


17.2 


-81 


932 


3.0 


6.8 


-3.5 


383 


S5 


0.0256 


0.04 


4000 


18.6 


-79.7 


876, 1147 


2.4 


6.9 


-3.2 


335 


S6 


0.0512 


0.04 


4000 


19.5 


-80.0 


815, 1064 


2.2 


6.9 


-3.2 


295 



Spontaneous firing with 6 components 

DRN SE and other raphe SE cells often fire spontaneously in a pacemaker fashion, 
so it is of interest to explore this aspect in a computational model, at first with just 
Ino,, Ikdr, In, Isk and Ibk and the leak current as in the last section. Because only 
small depolarizing currents were required to induce firing in the most of the runs 
of the previous section, it seemed that a small increase in excitability could induce 
spontaneous firing. This was achieved with just one change, namely an increase 
in the value of k^ ai from its prior value 8 to 10 mV. A set of results for fi = is 
given in Table 23. In all the runs of this section, which are given a P label, the V- 
dependent forms for the sodium and potassium time constants are used, the sodium 
conductance is gNa,max=0.Q75, gKDR,max=0- 0768 and gsK,max=0. 012. Further, here 
the activation parameters of 1^ are left the same, with = — 25mV and = 7. 

In the first of this set of runs, PI, the N-type Ca 2+ conductance is 0.012. The 
cell fires regularly with an ISI of 467 ms and the duration of the spikes is 2.3 ms. An 
increase in gN,max to 0.24, as in P2, increases the ISI to 982 ms and the maximum 
internal Ca 2+ to 512 nM. The effects of increases in cell capacitance C and area 
A are explored in P3 to P5. As a general rule, increasing the capacitance leads to 
increases in the ISI, spike duration and peak internal calcium concentration whereas 
increasing the area decreases the ISI and the peak internal calcium ion concentration. 

The simplified BK current model, as described above, is included in runs P6 
and P7. Compared to P2 which has the same C and A values, the ISI is shortened 
somewhat but most significantly the spike duration is reduced to 2.15 ms in P6 and 
2.0 in P7, these values being in the experimental range for DRN SE neurons. In 
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contrast to the runs S5 and S6 (induced firing), the spike trains in the pacemaker 
runs P6 and P7 with BK current were regular with a constant ISI. 



Table 23: Pacemaker firing with \x = and fc/v ai = 10, changes in C,A and effect of 
BK 



Run 


GN 


GBK 


C 


A 






ISI 


D 


INA 


IK 


Ca,i, max 


PI 


0.012 




0.04 


4000 


26 


-79 


467 


2.3 


8.4 


-4.0 


250 


P2 


0.024 




0.04 


4000 


27 


-84 


982 


2.45 


8.4 


-4.0 


512 


P3 


0.024 




0.05 


4000 


23.5 


-84.2 


1099 


2.8 


8.3 


-4.0 


554 


P4 


0.024 




0.04 


5000 


27 


-83.5 


790 


2.4 


8.5 


-4.1 


430 


P5 


0.024 




0.05 


5000 


23.3 


-83.8 


887 


2.75 


8.3 


-4.0 


445 


P6 


0.024 


0.0256 


0.04 


4000 


25.5 


-83.6 


830 


2.15 


8.4 


-3.3 


423 


P7 


0.024 


0.0512 


0.04 


4000 


23.9 


-83.0 


710 


2.0 


8.5 


-3.4 


362 



Low threshold calcium current, It 

The low threshold Ca 2+ current It has often been ascribed a role in pacemaking 
in DRN SE neurons (for example, Aghajanian and Sanders-Bush, 2002). In this 
subsection we investigate the effects of the inclusion of this component in the model 
for which some results without It were given in Table 23, making the number of 
components 7, including leak current. Experimentally this would correspond to the 
block of I a and In, the L-type Ca 2+ current, which is relatively small in these cells, 
being approximately combined with the N-type. 

Here the value of gT,max was estimated from voltage-clamp data in Penington et 
al. (1991), yielding a value of approximately 0.05 /iS for a cell at room temperature, 
this estimate therefore being considered low for a cell at body temperature. A value 
of gT,max=0.l25 //S is a reasonable upper limit, which compares favorably with an 
earlier estimate of 0.180 /xS from the voltage-clamp data in Burlhis and Aghajanian 
(1987) - but see also the next subsection where I a is included. 

When, in a run labeled Tl, the higher value (0.125) was employed for gT,max, 
using the same parameters as in P7 of Table 23 (with other conductances as in the 
runs (S) of Table 22), there was an early spike with a small notch at 15 ms, a doublet 
at 300 ms, and spike which did not repolarize for over 100 ms, with a resultant huge 
increase in internal Ca 2+ concentration. This contrasts to the regular train in P7 
with gT,max=Q- In run T2, with gT, ma x at the lower estimate of 0.05 /zS there were 
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Figure 17: Doublet spikes appear at regular intervals (about 1500 ms) when the It 
conductance is around 0.05 /iS in Run T2 of the 6-component model. 
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Figure 18: Consecutive spikes with an ISI of 734 ms when the It conductance is at 
0.025 fiS (Run T4). 

initially 4 short ISIs involving singlet spikes with concomitant Cai levels of about 
500 nM, followed by a regular train of doublets with an ISI of 1557 ms, V max = 28 
mV, then after the doublet V m i n = —84 mV, and a Cai rising to about 945 nM. The 
doublet spike is pictured in Figure 17 and has the appearance of doublets reported 
for DRN SE neurons by Hajos et al. (1996). 
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Because the internal Ca 2+ concentration became somewhat elevated, a stronger 
calcium pump (constant K s in Table 17) was tested, as in run T3, which made 
the spike train more irregular. Initially there were 6 singlets with notches of in- 
creasing magnitudes on their repolarizing edges, with V max values 22 mV, AHPs 
to -82 mV and durations about 2.1 ms. This was followed, at intervals of over 
2000 ms, in sequence, by a doublet, a triplet, and then several quartets which 
eventually degenerated into a damped oscillation. The internal Ca 2+ concentra- 
tion for the first spikes grew to about 400 nM, with later levels around 3400 
nM. Surprisingly, in Run T4, with gT,max reduced to 0.025, somewhat below the 
lower estimate, a very regular spike train occurred with the following properties: 
X = [25.5,-84,2,734,9.1,-3.5,462] and I T , max = 0.2 nA. The last three ISIs all 
had values of 734 ms. A pair of spikes is shown in Figure 18, indicating a small 
notch on the repolarization phase. 

With the same value of gT,max, the half-activation and half-inactivation poten- 
tials of It were increased by 5 mV in Run T5A and decreased by 5 mV in run T5B. 
In the former case, the spike train properties were similar to those described in the 
previous example, but the notch on repolarization was smaller. In the latter case, 
however, the membrane potential took on a new character between spikes as now 
there was a smooth broad undulation (notch) preceding the spike as V retreated and 
paused before swinging up to threshold, similar to the appearance of some spiking 
in DRN SE cells in Bayliss et al. (1997). There was no discernible notch on the 
repolarization phase, and the It current was noticeably smaller with a maximum 
value of only 0.03 nA. The ISI was about 702 ms. 

The next four runs, designated T7-T10, were designed to ascertain the relative 
roles of It and I^ a in spike triggering, especially in pacemaker mode. The higher 
(+5 mV) half- activation and half-inactivation potentials were employed for It- In 
T7, the smallest value of gT,max= 0.025 /zS was used with gNa,max smaller by 25% at 
0.54 fiS. Previously this value of the T-type conductance gave rise to bursting, but 
now the spike train was perfectly regular with an ISI of 642 ms and spike duration 
2.05 ms. Further, there was no sign of a notch either in the late part of the ISI or 
the repolarization phase, but the maximum of It was quite small at 0.028 nA. In 
T8, with the higher 5'T,maa;=0.05 and the same gNa,max, a regular spike train also 
ensued with an ISI of 705 ms, maximum It of 0.59 nA and with a small notch on the 
repolarization phase of the spikes. In run T9, as depicted in Figure 19, the sodium 
conductance was returned to its higher (by 25 %) value and with gT,max=0.05, there 
was now no bursting with an ISI of 802 ms, a small post-spike notch and a maximum 
It of 0.725 nA. Noteworthy is that the notch is higher at about -60 mV, compared 
with about -73 mV in the previous figure. It seems that the higher the notch, the 



66 



closer the cell is to bursting mode. 
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Figure 19: Spikes with an ISI of 802 ms in run T9 with higher half-activation and 
half-inactivation potentials. 

Finally, in run T10, the value of gNa,max was decreased to 0.405 /iS (0.75 x 0.54) 
and gr,max was increased to 0.075 /iS. In this case spiking halted after 5 spikes, the 
4th being a doublet and the 5th being a multiple oscillation accompanied by rises 
in Cdi to about 4000 nM. 

Inclusion of Ia- 8 component model 

A strong outward current, labeled as I a, was found in DRN SE neurons by Se- 
gal (1985) and studied in Aghajanian (1985), Burlhis and Aghajanian (1987) and 
Penington and Tuckwell (2012). Such a current has been postulated as exerting 
considerable control over spiking activity. Here we first describe two interesting re- 
sults obtained with fairly small Ia and It, and then obtain a very satisfying result 
by using conductances estimated from voltage-clamp data, in which I a increases in 
amplitude throughout the ISI and concomitantly It decreases, reinforcing the idea 
that competition between these two components is an important factor in deter- 
mining the ISI. Finally, in 11.8.3, a methodical analysis results in the culmination 
of computations of notch-free spikes which are often reported for these cells. 

Results With Small g A ,max, 9T,max 

In the first of these preliminary runs, called Al here, the following maximal con- 
ductances (in /iS) are employed: g N a,max=0.594:, g K DR,max=0.0768, g N)max =0.024:, 
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gSK,max=0 .012, g B K ,max=2 9 K D R,max / 3 , gA,max=gT,max=0.l. The activation and ill- 

activation parameters for I a and It are as given in Table 16, the value of the 
pump strength K s is fairly small at 0.3125 x 10~6 and Vn x is set at the rel- 
atively low value of -25 mV. The resulting voltage trajectory is shown in Fig- 
ure 20 and the regular train of spikes have the following characteristics. X = 
[16raV, —77, 2.1, 1185, 6.55, —2.75, 433], a small notch with magnitude about 1.5 mV 
at the bottom of the fast post-spike repolarization and no prolonged plateau at 
the end of the ISI. The other maximum current amplitudes are lT,max = 0.026, 

lA,min = —0.12 and In, max = 0.45. 




Figure 20: The membrane potential showing consecutive spikes with an ISI of 1185 
ms, with gA,max=gT,max=0.l, as described in 11.8.1. 

In the following run, called A2, only 1 change was made, namely to increase 
gN,max by 50% to 0.036, but this resulted in a radical change in the nature of 
the spikes, which are depicted in Figure 21. A very long plateau developed as 
the membrane potential lingered for over 2000 ms about 5 mV above resting level 
before an oscillation of growing amplitude called pre-spike oscillation, similar to 
that observed in some cells (Atherton and Bevan, 2005) finally triggered a spike. 
Apart from the very long ISI of about 3480 ms, the remaining spike properties such 
as maximum V, minimum V, spike duration and magnitudes of I^a, Ikdr, I a and 
It were practically the same as in the run with smaller gN,max- Figure 22 shows 
that during the long plateau, the internal Ca 2+ concentration remained practically 
constant at a low level, down from a somewhat higher peak value of 480 nM. 
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Figure 21: Spikes obtained when in the previous parameter set only one change was 
made, consisting in an increase in the N-type Ca 2+ conductance. 
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Figure 22: Time course of internal Ca 2+ concentration in the ISI during the spiking 
shown in the previous graphic. 

Further voltage clamp estimates of conductances for I a and It 

Using the updated model obtained with 7 components (plus leak), a careful re- 
analysis was performed for the voltage-clamp data in Burlhis and Aghajanian (1987), 
but using the linear form for the calcium current and not the constant field form. 
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Recall that two sets of clamps were performed, from -80 mV to -56 mV and from - 
65 mV to -56 mV. In this analysis, variations on the activation and inactivation 
parameters for It and I a, namely V^, Vr 3 , Va 1) Va 3 , were tested against the 
values given in Table 16, which were essentially those from Huguenard and Mc- 
Cormick (1992). There were 8 variations for each of the two holding potentials as 
each half-activation or half-inactivation potential was increased or decreased by 5 
mV. The conclusion was that the values in the Table could not be improved upon. 
The final values of the maximal conductances resulting from this analysis were (in 
/j,S) gA,max=0.75 and gT, ma x=0.265. These values were used in conjunction with 
the conductances g N a,max= 0.594, gKDR,max=0. 07 '68, g N , max =0.024:, g S K,max=0M2, 
gBK,max=2gKDR,max/^ and a small depolarizing current of /i = —0.1 nA. This set of 
parameters will be referred to as run A3. 

Results are shown in Figures 23 and 24. In the first of these Figures, the 
two top panels show the time course of the membrane potential V in mV and 
the internal Ca 2+ concentration Ca>i in mM versus time in ms. In the lower 4 
panels are plotted the currents I^ a , Irdr, and 1^ on expanded time scales and 
Isk all in nA. In the membrane potential trajectory a small notch (a few mV) 
can be seen at the base of the falling edge of the spike. The membrane potential 
climbs steadily to the threshold of the succeeding spike. The spike properties are 
X = [11.9,-74,2.2,1363,5.3,-2.3,300] and I N , max = 0.34nA I SK , m in = -0.18nA 
and lBK,min = —1-55 nA. In Figure 24 are shown the details of the currents I a and 
It during the ISI. According to these calculations, both currents peak abruptly dur- 
ing the spike but in the period until the subsequent spike I a gets progressively but 
slowly more negative, while It reaches a maximum at about 400 ms (in this exam- 
ple) post spike and then declines steadily until it jumps up to trigger the next spike. 
The competition between I a and It is therefore a factor in determining the length 
of the ISI, as postulated by many. Thus, these calculated activities can be com- 
pared with the proposed mechanisms of pacemaking in these cells as summarized 
by Jacobs and Azmitia (1992). 

During an action potential, calcium enters serotonergic neurons through a high- 
threshold calcium channel. This is followed by a large (15-20 mV) AHP generated 
by a calcium- activated potassium conductance. This AHP results in a long relative 
refractory period, thus preventing discharges in bursts and insuring slow rates of 
firing. As the AHP decreases (due to sequestration and/ or extrusion of calcium), it 
deinactivates a low-threshold calcium current and an early transient outward potas- 
sium current (I A). The currents generated by the activation of these two voltage- 
dependent channels are opposed, with IA tending to slow the rate of depolarization 



70 



and the low-threshold calcium conductance increasing it. Under normal conditions, 
the calcium conductance is stronger, thus leading to a shallow ramp depolarization, 
which ultimately reaches threshold, fires, and, as the calcium enters the cell, reini- 
tiates the sequence of events. The slope of this ramp is what determines the rate of 
discharge of serotonergic neurons. 
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Figure 23: Time courses of several variables versus time during spiking and during 
the ISI for Run A3. Top left, membrane potential V, top right, internal Ca 2+ 
concentration with a maximum of 300 nM. In the bottom four panels, from left 
to right, are shown the fast sodium current, the high threshold calcium current, 
the delayed rectifier potassium current and the calcium-dependent (SK) potassium 
current. Note the different time scales for the faster and slower currents. 



Notch-free spikes 

As seen above, many spikes generated by the model, including those in the last 
subsection, which have taken account of the major known currents operating in 
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Figure 24: The computed time courses of It (top panel) and I a (bottom panel) 
during and beween spikes in the model with 7 component currents plus leak for 
Run A3. 

DRN SE cells, have notches of various magnitudes, appearing usually on the falling 
edge of the spike. Such notches are likely to have their origin in equilibrium points 
such as spiral points in the 8 or 9-dimensional space in which the trajectories move. 
Notches in DRN SE spikes are not uncommon - see for example Hajos et al. (1996). 
The majority of reported spikes, however, are notch free, or possibly with such small 
notches they are not noticeable. 

As a first step towards finding the factors which might lead to a notch or its 
absence, the parameters and Vr 3 ,, being the half-activation and half-inactivation 
potentials for It, were varied over several mV from their standard values. The notch 
size with the standard values ws 4.5 mV at the base of the spike, as seen in Figure 
23. All the tested values resulted in similar sized notches, so it was concluded that 
notch size was practically independent of the magnitudes of these two parameters. 

Secondly the parameter gsK,max was varied from its standard value of 0.012. 
Smaller values of gsK,max gave rise to larger notches, up to 7 or more mV, whereas 
larger values gave smaller notches and also in some cases an undulating plateau. 
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Thirdly, the strength of the calcium pump, through the parameter K s in Table 
17, was varied. Increases in this parameter led to larger notch sizes and faster 
spiking. With smaller values, notch size also grew to over 8 mV and spiking became 
extremely slow. 

Finally, the maximal conductance gKDR,max was varied from its standard value 
of 3 x 0.0256. The notch sizes for gKDR,max = 2-5 x 0.0256, 2.0 x 0.0256 and 1.75 x 
0.0256 were respectively 4.5, 2.7 and 0.15 mV indicating that the notch size varied in 
an inverse fashion to the delayed rectifier conductance. Somewhat amazingly when 
the value of gKDR,max was set at 1.5 x 0.0256, the spikes subsequent to the first, 
which had a very small notch of about 0.15 mV, became absolutely notch free. The 
results of such a run, called A4, are depicted in Figure 25, the quantities appearing 
therein being the same as in Figure 23 for run A3. 
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Figure 25: Voltage, current and internal Ca 2+ concentration (as in Figure 23) for the 
regular spike train in Run A4 where notch-free spikes were obtained after reduction 
of the conductance of the potassium delayed rectifier current. 



73 



12 Spontaneous activity with and without J# 



The question as to whether DRN SE neurons require excitatory synaptic input to 
fire regularly in pacemaker fashion has been much discussed. It is generally expected 
that in slice the amount of synaptic input would be small or zero, so the fact that 
there are spontaneously active neurons in slice indicates that synaptic input may 
not be required. Note, however, that Segal (1985) shows a neuron in slice with 
definite EPSPs. 

In this section we focus on modeling the possibility of regular spiking with no 
external drive so that /i — 0. It was not so easy to achieve this without I H , the 
general finding being that small changes in many parameters would frustrate the 
occurrence of regular firing, leading to such results as bursting either in doublets or 
higher multiplets and long duration spikes with concomitant enormous increases in 
internal Ca 2+ concentration. 

There are doubtless many combinations of parameters which give rise to regular 
spontaneous actuivity without the introduction of Ih and here we report on one such 
set. The significant features of this set are (a) there is no Ih nor any 1^ (b) the fast 
sodium conductance is relatively much larger at gNa,max=l-35 and VVai = —33.1 
mV is at its standard value (c) Vjvi = —25 mV (d) the magnitudes of gA,max=®-3 
and gT,max=0.15 are considerably reduced (e) the area and capacitance have their 
standard values of A=4000 and C=0.04 (f) the calcium pump strength is K s = 
0.625 x 10~6 (g) gKDR,max=3 x 0.0256 (h) the parameter CSF which determines the 
effective increase in Cdi due to Ca 2+ influx from 1^ is decreased to 0.25 compared 
with the standard value of 0.7. 

With these changes spontaneous activity ensued but did not settle to regular 
firing until after over 20 spikes. When regular (fast) spiking was established, with 
an ISI of only 282 ms, the spikes had the following properties 

X = [34, -77, 1.9, 282, 16, -4.05, 248]. 

Further, a large post-spike notch with an increment of 8.5 mV occurred at -71 mV. 

Spontaneous activity in the complete simplified model, including I L 

The presence of Ih makes it easier to achieve spontaneous spiking activity. First we 
note that diverse spike patterns can be obtained with various calcium pump rates for 
certain parameter sets, as illustrated in Figure 26. Here the values of the 
pacitance and maximal sodium conductance are set at higher than standard values, 
being 0.08 nF, 5200 sq microns and 1.188 /iS. Remaining conductances, in /xS are 
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9KDR,max = 1 • 5 x 0.0256, # T ,max=0.95 x 0.265, g A , ma x=0.9 x 0.75 and g IH}max =0M2. 
Also, in mV, V Na ,i = -36.1, k NaA = 10, V Nat3 = -53.3, VT, 1 = -57, VN1 = -25 
with other parameters taking standard values. 

When the pump strength is K s = 0.25 x 1.25 x 10~ 6 , the spike train is a regular 
sequence of singlets with an ISI = 2102 ms and duration 3.2 ms, as depicted in the 
top record of Figure 26. The level of Cdj reaches 575 nM. A somewhat larger pump 
strength of K s = 0.35 x 1.25 x 10 -6 leads eventually, as seen in the middle panel of 
Figure 26, to a regular pattern of spikes with an ISI of 1501 ms, but initially there is 
a short interval terminating with a doublet and then a much longer interval before 
the train settles to its periodic form. During the latter part of the ISI a small kink 
is seen in the voltage trajectory. Similar kinks have been observed in spike trains 
of some DRN SE neurons (Park, 1987; Li and Bayliss, 1998). 

At a still higher pump strength K s = 0.40 x 1.25 x 10 -6 , the spike train becomes 
more irregular, settling to alternating long and short ISIs of durations about 1500 
and 1000 ms, respectively. The short ISIs display no kink but the longer ones have 
a pronounced kink about 75% into the ISI, similar to those found in some midbrain 
dopamine neurons (Neuhoff et al., 2002). Initially there is a very short interval 
which in this case terminates with a triplet. 

50 



-50 
-100 

1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 

50 


V(mV) 

-50 
-100 

1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 
50 1 1 1 1 1 1 1 1 1 1 



-50 

-100' ' 1 1 1 ' 1 1 1 1 1 

1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 

t(ms) 

Figure 26: Spike trains in the model with Ih, but not II, with three different 
calcium pump strengths, K s , in multiples of 1.25 x 10~ 6 . Note the presence of a 
doublet (D) and a triplet (T) with the stronger pumps. 
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Spontaneous firing with only small changes in parameters from those for driven 
activity 

The properties of the spikes in runs A3 and A4 were close to those of DRN SE 
neurons, but in those runs there was an applied depolarizing current of magnitude 
\x = —0.1 nA required for spiking. It was found that with very minor changes to, 
for example, only 5 parameters together with the introduction of In, spontaneous 
pacemaker-like activity ensued. The very small changes made were, as summarized 
in Table 24, 5% in the half- activation potentials for I Na , I T and I a, as well as a 
5% change in the slope factor for the steady state sodium activation. Further, the 
maximal sodium conductance, gNa,max, was reduced to 0.955 of its value. There are 
three sets of runs here based on this slightly altered set, labeled SS, SSL and F. 



Table 24: Parameters which differ for spontaneous and driven activity 



Parameter 


Driven (A3) 


Spontaneous (SSI) 




-33.1 


-34.755 




10 


10.5 


v Tl 


-57 


-54.15 




-60 


-57 


9Na,max 


0.594 


0.567 


9 1 H, max 





0.012 


fi 


-0.1 






Table 25: ISIs of pacemaker spikes with I H but not II 



Run 


9lH,max 


v Nl 


Vn 3 


K s 




9N,max 


ISI 


SSI 


0.012 


-25 


-50 


0.25 x 1.25 x 1(T 


-6 


0.024 


1950 


SS2 





-25 


-50 


0.25 x 1.25 x 10~ 


-6 


0.024 


No spikes 


SS3 


0.012 


-10 


-35 


0.25 x 1.25 x 10- 


-6 


0.024 


1200 


SS4 


0.012 


-10 


-45 


0.25 x 1.25 x 10- 


-6 


0.0462 


1640 


SS5 


0.012 


-10 


-45 


1.25 x 1.25 x 10- 


-6 


0.0462 


1500 



Before describing results with all components, we consider the set (SS) in which 
pacemaker activity is observed with Ih but II is not included. There are 5 such 
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runs described in Table 25. Table 24 gave the values of the parameters which 
differ between SSI (spontaneous) and A3 (driven). In Table 25 there are only 5 
parameters, giH,max, Vn 1 , Vn 3 , K s and gN,max, which are varied. The half-activation 
potential of In could be higher than -25 mV, so the value -10 mV was tested. 
However, in either case In is only activated substantially above action potential 
threshold so this change probably only has a minor effect. The value of gN,max was 
also estimated from voltage-clamp data in Penington et al. (1991) to give the value 
0.0462 fiS, which is the value used in runs SS4 and SS5, but again, this change 
had only small effects. The actual value in native cells at body temperature is not 
available. 

For these parameter sets, the computed model ISIs ranged from 1200 to 1950 
ms. All the spike trains were regular, and only in the last run of this group, SS5, was 
there a small kink about half way along the ISI (cf Park, 1987; Li and Bayliss, 1998; 
see also spiking in mouse locus coeruleus, de Oliveira et al., 2010). All other spike 
trains had properties resembling those often observed in rat DRN SE neurons, which 
we note may have frequencies of discharge as low as 0.1 to 0.25 Hz (Aghajanian et 
al., 1978; Hajos et al.,1995 ; Bambico et al., 2009) corresponding to ISIs as large as 
2500 to 10000 ms. 

Spontaneous firing in the complete simplified model 
The complete simplified model has the 10 currents 

I A, IkDR, It, II, In, Ih, Inu, IsK, IbK, Iheak, 

the M-type potassium current being omitted here, but to be included in future work. 
It is rather worth pondering that such a complex system, ofl6 nonlinear differential 
equations, admits smooth and regular periodic solutions, especially in light of the 
fact the biological existence of such periodic solutions is essential for the proper 
functioning of many neurons in the mammalian and other nervous systems. In the 
calculations described so far, it was apparent that the ranges are narrow of all of 
the 90 parameters in the computational model in which regular periodic spiking 
occurs. For the proper functioning of the neurons and hence the brain, this implies 
there is much scope for the occurrence of serious problems in maintaining required 
or normal activity. 

Altogether about 30 different parameter sets in the 10-component current model 
were used to explore the roles of some parameteras, of which 24 will be described 
briefly. These are the two sets SSL (i.e., SS with added L-type) and F (for final) 
there being 14 sets in the SSL-sequence and 10 in the F-sequence. In the SSL set, 
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for which results for some spike properties are given in Table 26, the surface area 
of the cell was always taken to be the standard value 4000 sq microns, except in 
the last two runs SSL13 and SSL14, and giH,max was always 0.012 fiS. There are 
only seven parameters which vary in this set. Thus an idea, albeit limited, of how 
the spike train is affected by changes in the parameters of the two higher threshold 
Ca 2+ currents, In and II as well as the Ca 2+ pump rate and area can be obtained. 

Generally larger pump rates make the ISI shorter (compare SSL2 with SSL3 
and SSL6 with SSL7 and SSL8) and larger values of gN,max make the ISI longer 
(compare SSL4 with SSL5). This is because higher pump rates clear internal Ca 2+ 
faster which has less chance to activate Isk channels and higher gN,max (or gL,max) 
leads to more calcium influx with subsequent activation of Isk- Increasing the half- 
activation potentials of II and In made the ISI smaller (compare SSL6 with SSL9), 
because there was less calcium influx to activate Isk- Increasing the calcium pump 
strength can induce repetitive firing which is absent for smaller pump strengths 
(compare SSL10 and SSL11). Comparing SSL12 and SS13, an increase in area A 
causes a speeding up of the spiking, but in comparing SSL13 and SSL14, there is a 
trade-off between the inhibitory effect of a weaker calcium pump and the excitatory 
effect of a larger area, the former dominating so that the ISI increases. 

The final set of results, labeled F, is given in Table 27. Here there are only 
4 parameters which vary, the following parameters being fixed: gL,max=0- 00462, 
gN,max=9gL,max, g i h ,max=0 .012 , and the half-activation and inactivation potentials 
for I L and I N , as in the last row of Table 26. The parameters that varied were 
the area, A, the pump strength, K s , and the sodium and It conductances gNa,max 
and gr,max- ISIs ranged from a smallest value of 506 ms to a largest value of 1694 
ms. A few spike trains (F2, F3) were irregular but mostly they were periodic. 
In the ISI there were instances of kinks, changes of slope (mainly not great) and 
plateaux, but no notches. The spikes in runs F4 ,F6, F7, F9 and F10 were smooth 
as usually observed for DRN SE neurons, and the spike trains were regular periodic. 
Investigations in which gx,max varied indicated that this component was probably 
the likely origin of kinks when they occurred. 

From the results in Table 27, comparing Fl and F2 again demonstrates the 
slowing down of the spike train by a decreased calcium pump strength. Decreasing 
the area A, as from F2 to F3, can change the character of the spike train from 
steady state regular to an alternating long-short sequence. Comparing F3 with 
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Table 26: ISIs of spontaneous pacemaker spiking in the complete model: SSL sequence 



Run 




Vl 3 




v Nl 


Vn 3 


Q TV TYICIX 


K s 


ISI 


Remarks 


SSL1 


-15 


-25 


0.01 


-25 


-50 


0.024 


3.125 x 10~ 7 


3000 


Plateau 1000 ms 

Ca itTnax 510 nM 


SSL2 


-20 


-50 


0.01 


-10 


-45 


0.0462 


3.125 x 10~ 7 


2757 


Similar to SSL1 


SSL3 


-20 


-50 


0.01 


-10 


-45 


0.0462 


3.90625 x 10~ 7 


1819 


Slope change mid ISI 
Stronger Ca pump 


SSL4 


-20 


-50 


0.01 


-10 


-45 


0.0693 


4.6875 x 10~ 7 


2700 


Late rise in V 
9N,max 50% larger 
Ca i)max 1100 nM 


SSL5 


-20 


-50 


0.01 


-10 


-45 


0.0231 


4.6875 x 10~ 7 


1114 


Kink in V 

Qn max 50% smaller 

Ca^ max 380 nM 


SSL6 


-20 


-50 


0.01 


-10 


-35 


0.024 


3.125 x 10~ 7 


1980 


Smooth trajectory 
dV/dt steadily decreases 
during the ISI 


SSL7 


-20 


-50 


0.01 


-10 


-35 


0.024 


4.6875 x 10~ 7 


1200 


Kink mid ISI 

QsKmax increased 25% 


SSL8 


-20 


-50 


0.01 


-10 


-35 


0.024 


3.125 x 10~ 7 


1932 


Slope change mid ISI 
9sK,max as in SSL7 


SSL9 


-15 


-50 


0.01 


-5 


-35 


0.024 


3.125 x 10~ 7 


1476 


Slope change mid ISI 
9sK,max standard 

and Vjvi 5 mV higher 


SSL10 


-25 


-50 


0.01 


-10 


-45 


0.0362 


3.125 x 10~ 7 




1 spike only 
Larger g N , max 


SSL11 


-25 


-50 


0.01 


-10 


-45 


0.0362 


3.90625 x 10- 7 


2400 


Stronger pump 


SSL12 


-20 


-45 


0.00462 


-10 


-45 


0.04158 


3.90625 x 10~ 7 


1659 


9 L, max- 9 N, max 1-9 


SSL13 


-20 


-45 


0.00462 


-10 


-45 


0.04158 


3.90625 x 10~ 7 


1346 


Slope change mid ISI 
A=5000 /i 2 


SSL14 


-20 


-45 


0.00462 


-10 


-45 


0.04158 


3.125 x 10~ 7 


1454 


Slope change mid ISI 
A=6000 fx 2 



F4, a decrease in calcium pump strength can change an irregular spike train into a 
regular one and a small change in this parameter can introduce a kink in the voltage 
trajectory during the ISI. A comparison of Fl and F6 shows that a small increase 
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in gNa,max and a concomitant decrease in gr,max can produce smooth trajectories 
between spikes rather than ones with a marked change of slope, but that the ISI 
is hardly changed. From F6 to F7 a decrease in area is inhibitory as well as the 
reduction in gT,max, but the spike train remains smooth and regular. F8 and F9 
demonstrate the speeding up of spiking with increased area, and in F9 an increase 
in gNa,max and a concomitant decrease in gr,max leads to a plateau towards the end 
of the ISI. Comparing F9 and F10 further demonstrates how a larger area leads 
to a smaller ISI. Finally, in FURA-2AM measurements under various conditions, it 
was found that the resting internal Ca 2+ concentration was about 32 nM, but when 
this value was employed rather than 50 nM, the spike properties differed very little, 
including ISI and the maximal Ca,. 

A set of computed results for a run very similar to F7 is given in Figures 27 and 
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Figure 27: Results showing spontaneous pacemaker activity (/x = 0) for the complete 
model in a run very similar to F7 in Table 27. All parameters are given in Table 
28. Top two panels show the membrane potential and internal Ca 2+ concentration 
during an ISI. The bottom left 3 panels show the details of the three currents Ij^ a , 
In, Irdr during a spike and the remaining panel shows the calcium-dependent 
current Isk during the ISI. 
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Table 27: Results for ISIs for pacemaker spiking in the complete model, F sequence 



Run 


A 


K s 






QNa,max 


9T,max 


ISI 


Remarks 


Fl 


6000 


3.90625 x 


10" 


-7 


0.567 


0.265 


1148 


Slope change mid ISI 


F2 


6000 


6.25 x 10- 


-7 




0.567 


0.265 


506 


Irregular start 


F3 


4000 


6.25 x 10- 


-7 




0.567 


0.265 


500,1200 


Alternate short, long 


F4 


4000 


5 x 10~ 7 






0.567 


0.265 


1500 


Regular train 


F5 


4000 


5.46875 x 


10" 


-7 


0.567 


0.265 


1300 


Kink 


F6 


6000 


3.90625 x 


10" 


-7 


0.594 


0.22525 


1145 


Regular smooth spikes 


F7 


4000 


3.90625 x 


10" 


-7 


0.594 


0.1855 


1694 


Regular smooth spikes 


















Ca iimax 550 nM 


F8 


6000 


6.25 x 10- 


-7 




0.594 


0.22525 


719 


Slope change mid ISI 


F9 


6000 


6.25 x 10- 


-7 




0.675 


0.14575 


770 


Plateau at end of ISI 


F10 


4000 


6.25 x 10- 


-7 




0.675 


0.14575 


1157 


Small slope change end ISI 
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Figure 28: Results for the same parameter set as in Figure 27. The currents It, I a 
and Ih are shown during an ISI and the currents II and Ibk are shown during a 
spike. 
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28. The currents which remain extremely small for most of the ISI, namely I^ a , 
Irdr, II, In and I B k are shown only during the spike, whereas the remaining 4 
currents are shown throughout the ISI. Again, the interplay between It and I a is 
seen to be important in that both are activated in similar voltage ranges and tend 
to oppose one another. No measurements have been made of any of these currents, 
during either a spike nor during the ISI, with which to compare the numerical 
results. 

Summary of parameters for spontaneous activity 

We give here a complete list of the parameters for the run F7 of Table 27 for which 
spikes were regular with an ISI of 1694 ms corresponding to a frequency of 0.59 Hz, 
as depicted in Figures 27 and 28. This parameter set, based as much as possible 
on experimental data, may be used as a starting point for more complex models 
with a more complete knowledge of the characteristics of each component, though of 
course there is always large amount of variability from cell to cell and under different 
conditions so there are no definitive sets of parameters. In Table 28, corresponding 
to the initial Table 16, are given the parameters for activation (and inactivation if 
appropriate) of the 9 channel types I Na , Irdr, It, h, In, I a, Ih, Isk and I BK . 
Table 28 differs slightly in format from Table 16, as parameters for the simple I BK 
current have been added and the time constants for sodium are V-dependent. Table 
17 of cell properties, calcium dynamics and equilibrium potentials is the same for 
F7, except for the one parameter for the calcium pump strength, which is now 
K s = 3.90625 x 10~ 7 . Table 29 gives the maximal conductances, all of which differ 
from those in Table 18, and there is the additional conductance for I BK . In all, 
there are 90 parameters as Faraday's constant is not a parameter. 
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Table 28: Activation and inactivation parameters for run F7 
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Table 29: Maximal conductances 



9Na,max 


0.594 


9N,max 


0.04158 


9KDR,max 


0.0384 


9A,max 


0.75 


9T,max 


0.22525 


9lH,max 


0.018 


9L,max 


0.00462 


9SK,max 


0.012 


9BK,max 


0.0256 







13 Anodal break 

Some experimental results on DRN SE neurons are able to be compared with model 
predictions. In Burlhis and Aghajanian (1987), some anodal break experiments were 
reported in which a sustained hyperpolarizing input current was suddenly removed, 
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with and without TTX. 
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Figure 29: Computed model anodal break voltage response with (assumed zero fast 
sodium conductance) and without TTX to a hyperpolarizing current applied for 
100 ms. 

For the former case, with gNa,max=0, in simulations a hyperpolarizing current 
fi = 0.1 nA was applied for 100 ms using the same parameter set as in run F10 
of Table 27. The result is shown in the left panel of Figure 29 where it is seen 
that during the application of the hyperpolarizing current, V decreases from resting 
level of -60 mV to about -79 mV with a relatively small time constant of about 7 
ms, which is similar to the average value reported in Li et al. (2001). In another 
run a time constant of about 30 ms was obtained, which is more in the range of 
reported values for the majority of DRN SE neurons. When released from the 
hyperpolarizing current, the voltage did swing up past rest to about -56 mV and 
eventually returned to resting level, mimicking approximately the result in Figure 
3A of Burlhis and Aghajanian (1987). The same experiment without TTX produced 
a spike on release as depicted in the right hand panel of Figure 29, with a large notch 
at about -63 mV. 
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14 f/I curves 



There do not appear to be published explicit frequency versus magnitude of applied 
current in DRN SE neurons. There are data on firing rate versus concentration of 
norephinephrine and phenylephrine (Judge and Gartside, 2006). For example, firing 
rate increased from about 1.25 Hz to 1.6 Hz when 10 /i M phenylephrine was applied. 
The computational model developed here was used to find frequency versus applied 
current (f/I) curves for two sets of parameters. In one set the cell fired spontaneously 
and even with a very small hyperpolarizing currrent, and in the other case the cell 
only fired if a depolarizing current of at least 77 pA was applied. The computed 
f/I curves are shown in Figure 30. In both cases, the firing rate increases as —\l 
increases as expected; but it was not expected that the firing rate would decrease 
slightly at values of the depolarizing curren 0.45 nA in the spontaneously active cell 
and 0.55 nA in the non-spontaneously active cell. Furthermore, just above these 
values the behaviors were quite different as the spontaneously active cell entered a 
regime of very rapid firing (76.9 Hz) of single spikes with fi = —0.5 nA whereas the 
non spontaneously active model neuron at /x = —0.6 nA began to fire in rapid but 
regular triplets. 

1.4 
1.2 
1 

f, Hz 

0.8 
0.6 
0.4 
0.2 



-6.1 0.1 0.2 0.3 0.4 0.5 0.6 

I = nA 

Figure 30: Calculated f/I curves for the cases of a spontaneously active neuron and 
a neuron with a threshold for firing of 77 pA. Note the slight downturns at large 
depolarizing currents. 
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15 Discussion 



Serotonergic neurons of the dorsal (and other) raphe have complex electrophysiology 
and neurochemistry about which much is known (Jacobs and Fornal, 1995; Azmitia 
and Whit aker- Azmitia, 1995; Aghajanian and Sanders-Bush, 2002; Harsing, 2006; 
Lanfumey et al., 2008; Lowry et al., 2008; Hale and Lowry, 2011) including their 
role in the therapeutical effects of SSRIs and other antidepressants (Piheyro and 
Blier, 1999; Guiard et al., 2011; Haenisch and Bonisch, 2011). 

Nevertheless, it is a remaining challenge to determine the nature and function- 
ing of the networks in which these neurons are involved, including their inputs from 
both nervous and endocrine systems and the effects they have on the many cells they 
influence. Considering the important roles played by these neurons in the physiolog- 
ical and psychological functioning of most mammalian species, it is not surprising 
that a more complete understanding is being pursued in various laboratories. Many 
of the details of the projections of DRN neurons and their interconnections have re- 
cently been determined (Vasudeva et al., 2011; Waselus et al., 2011), some with the 
use of genetic techniques (Bang et al., 2012). It is expected that quantitative models 
of the systems affected by and affecting DRN SE neurons will play an important 
role in confirming hypotheses concerning their involvement. 

Opposing views have been stated concerning whether DRN SE cells' apparent 
pacemaker activity is autonomous or does in fact require some external excitatory 
influence. There are many central neurons that have apparent spontaneous pace- 
maker activity, but in each cell type the mechanisms are seemingly different. In 
midbrain dopamine (DA) neurons, particularly of the substantia nigra in rats (Kang 
and Kitai, 1993; Wilson and Callaway, 2000; Durante et al., 2004; Atherton and 
Bevan, 2005; Foehring et al., 2009) and mice (Puopolo et al., 2007), L-type calcium 
currents, activated at relatively low membrane potentials have been found to be pri- 
marily involved though other voltage-gated calcium currents are present (Cardozo 
and Bean, 1995). The oscillations are set up primarily as an interaction between the 
L-type Ca 2+ and the calcium-dependent potassium currrent, with significant effects 
on frequency and regularity of spiking by SK3 channel activity (Wolfart et al., 2001), 
I A (Liss et al., 2001), buffering (Neuhoff et al., 2002; Foering et al., 2009) and I H 
(Neuhoff et al., 2002). The latter current is also important in pacemaking activity in 
many neuron types (Pape, 1996) including some hippocampal and thalamo-cortical 
neurons (Maccaferri and McBain, 1996; Luthi and McCormick, 1999). 

Locus coeruleus neurons seem to have quite different pacemaking mechanisms 
from those in DA neurons. In rats, these cells are spontaneously active both in vivo 
and in slice with frequencies between 0.4 and 4 Hz (Aghajanian et al., 1983; Alreja 
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and Aghajanian, 1991). Calcium-dependent potassium currents are the principal 
factor in the post-spike afterhyperpolarization and frequency of firing (Aghajanian 
et al., 1983) and availabilty of cAMP is required for sustained firing (Alreja and 
Aghajanian, 1991). In mice, de Oliveira et al. (2010, 2011) found that pacemaking 
in the majority of locus coeruleus cells was mainly due to sodium and potassium 
currents and that calcium currents were negligible during the ISI. There are differ- 
ences between rat and mouse, and development is also a factor (de Oliveira et al., 
2011). 

Other examples of cells exhibiting pacemaker-like activity include thalamic neu- 
rons (Mulle et al., 1986), cells in the suprachiasmatic nucleus (Jackson et al., 2004; 
Brown and Piggins, 2007; Itri et al., 2010), brainstem neurons involved in the con- 
trol of respiration ( Richerson, 1995; Rybak et al., 2003; Ptak et al., 2009; Dunmyre 
et al., 2011; Ramirez et al., 2011) cardiac cells (Brette et al., 2006; Mangoni et al., 
2003; Benitah et al., 2010; Vandael et al., 2010) ) and uterine cells (Rihana et al., 
2009). Further examples are summarized in Ramirez et al. (2011), which contains 
a detailed discussion of pacemaking in preBotzinger complex neurons where Ih and 
a persistent sodium current, Ino,p, play significant roles. 

Turning attention to DRN SE neurons, it has long been recognized that they 
undergo slow and regular spiking in a variety of conditions in most mammalian 
species with a rate of about 0.5 to 2 per second in the quiet waking state, higher in 
waking arousal, less in slow-wave sleep and zero in REM sleep (Jacobs and Fornal, 
1995; Aghajanian and Sanders-Bush, 2002). Such rates may also be broadly species 
dependent. The result of such regular firing is a fairly constant rate of release of 
serotonin in the many target areas (Jacobs and Fornal, 1995), both local and remote, 
but that the rate depends on physiological state. Since, according to Azmitia and 
Whitaker-Azmitia (1995), nearly every neuron in the brain is close to a serotonergic 
fiber and influenced either synaptically or by volume diffusion by serotonin, the 
consequence of the state-dependent release is a modulation of the activity of the 
whole brain activity as well as effects at spinal targets such as motoneurons. 

During regular almost periodic firing, there is observed the sequence of current 
flows already alluded to in previous sections and satisfactorily predicted by the 
computational model. Such a sequence has been described by many authors (e.g., 
Burlhis and Aghajanian, 1987; Jacobs and Azmitia, 1992; Piheyro and Blier, 1999). 
After a spike there is hyperpolarization due to calcium-dependent potassium cur- 
rent, followed by a fairly slow ramp-like approach to threshold, often with a plateau, 
as internal calcium is cleared and I a and It compete, with It assisting in or leading 
to the threshold crossing as I a retards the approach. 

Concerning whether DRN SE neurons are endogenous pacemakers whose regular 
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firing or pacemaker-like activity does not require any external stimulation, some 
early reports stated that there were only minor differences between in vitro and 
in vivo firing rates (Mosko and Jacobs, 1976; Jacobs and Azmitia, 1992) which 
supported the claim of intrinsic mechanisms for pacemaking. Some authors were 
skeptical (Crunelli et al., 1983) whereas many expressed the view that the ramp- 
like approach to threshold was sufficient evidence of the non-necessity of external 
influences (Jacobs and Azmitia, 1992; Aghajanian and Sanders-Bush, 2002). 

However, Vandermaelen and Aghajanian (1983) had found that most cells in 
vitro were less active than in vivo, which was supposed due to reduced excitatory 
drive. Pan et al. (1990) and Kirby et al. (2003) found only a small fraction of cells 
were active in slice, the latter authors claiming that such observations were probably 
evidence that noradrenergic input was cut off (see also Ohliger et al., 2003). 

Convincing physiological evidence for a role of noradrenergic drive was presented 
by Baraban and Aghajanian (1980) who also provided anatomical evidence that no- 
radrenergic input directly synapsed on DRN SE neurons (Baraban and Aghajanian, 
1981). The physiological evidence was that a-adrenoreceptor antagonists reduced 
the firing rates, although there was a possibility that the DRN SE neurons were sub- 
sequently able to recover and revert to their usual activity. Further, locus coeruleus 
cells are silent during REM sleep, during which DRN SE cells are also silent. Levine 
and Jacobs (1992) also noted that there was a dense noradrenergic input to the DRN 
and that the lack of effect of iontophoretic applications of noradrenergic agonists 
was possibly due to the already maximal noradrenergic input from locus coeruleus. 
Nevertheless, noradrenergic antagonists did suppress firing in anesthetized animals, 
but interpretation of this effect is complicated by the effects of some of these drugs 
on sodium channels. 

The noradrenergic excitation of DRN SE cells was established as being via by 
ax-adrenoreceptors and the mechanism of action of ax-adrenoreceptor agonists was 
explored by Aghajanian (1985) who found that I a was reduced by them, resulting 
in shorter plateaux and faster firing. However, with noradrenergic inputs there 
are uncertainties, as opposing effects can be elicited by the stimulation of different 
receptor subtypes, such as the ct\- and a2-adrenoreceptors (Pudovkina et al., 2003; 
Bortolozzi et al. (2003); Pudovkina and Westerink, 2005; Harsing, 2006; O'Leary 
et al., 2007). 

Since in brain slice preparations, most of the afferent inputs are removed, then 
the probable loss of most of the noradrenergic excitation is only one factor leading to 
a change in firing rate. As one example, inhibitory inputs, such as those originating 
in forebrain but mediated locally, will also be cut off, which would compensate to 
an unknown degree the loss of noradrenergic excitation. 
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In vitro preparations are also likely to have reduced or removed excitatory inputs 
from orexinergic neurons in the lateral hypothalamus (see Ohno et al., 2008 for 
review) and from histamine neurons as well as noradrenergic inputs which would 
contribute to their silence (Brown et al. 2002). Orexin afferents, which are dense 
in the DRN, affect both serotonergic neurons and GABA-ergic neurons such that 
at low levels the effect is primarily excitatory but at higher levels is inhibitory (Liu 
et al., 2002). The effects of orexinergic input on DRN SE cells are dependent on 
the state of arousal of the animal (Takahashi et al., 2005) having little effect during 
wakefulness but significant effect during REM and slow wave sleep. There are also 
differences in the strength of orexinergic inputs received by the dorsal raphe and 
medial raphe nuclei (Tao et al., 2006), with orexin A strongly targeting the DRN. An 
interesting feedback loop exists between the orexinergic and serotonergic systems, 
as the latter inhibit orexin neurons through 5-HTi^ receptors and associated GIRK 
currents (Ohno et al., 2008). Del Cid-Pellitero and Garzon (2011) have reported that 
there are loops involving the mPFC, the DRN and hypothalamic orexin neurons, 
which further complicates the matter of in vitro versus in vivo firing of DRN SE 
cells. 

Also relevant is the dense innervation of the DRN by CRF fibers as mentioned 
in the Introduction. Experiments using either urocortin or CRF have shown that 
at small doses SE cells are inhibited but at high doses excited (Kirby et al., 2000; 
Ordway et al., 2002; Pernar et al., 2004; Abrams et al., 2004). However Lowry et al. 
(2000) found that in vitro, application of CRF led to rapidly increased firing in DRN 
SE cells, but only in the ventral DRN. Additionally, CRF excites locus coeruleus 
neurons (Ordway et al., 2002; Snyder et al., 2012) which in turn excite many DRN 
SE neurons. Another relevant effect on DRN SE cells is from substance P neurons 
which excite them through a high density of receptors in the DRN (Ordway et al., 
2002). In conclusion, there are many factors involved in the differences between in 
vivo activity and in vitro activity so that a generalization is not possible. The most 
important factor is probably the variation of input patterns from one location to 
the other (Abrams et al., 2004; Lowry et al., 2008). One source of variability in the 
findings in different experiments could be the concentrations of ions such as K + , 
Mg 2+ and Ca 2+ in the extracellular fluid. 

The computations reported in the present article constitute a first step in the 
development of a quantitative description of the firing activity of DRN SE neurons. 
It has been demonstrated that these neurons may fire in a regular pacemaker fashion 
with or without any external drive and that only small differences in the properties 
of some ion channels can lead to one or the other mode of pacemaker activity. The 
present results are preliminary because of uncertainties in the properties of many 
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of the component currents. Further, simplifications and omissions have been made 
in the model in order to obtain some insight into the behavior of these neurons 
with their very complex set of ion channels and inputs. Apart from more elaborate 
descriptions of some of the dynamics of some of the component currents, future work 
will involve additional exploration of the changes in spiking activity as functions of 
changes in many of the parameters as well as including the effects of the various 
neurotransmitters affecting the activity of DRN SE neurons, as alluded to in the 
present discussion. One important transmitter is serotonin itself, especially through 
its activation of autoreceptors that reduce the firing rate. However, the serotonin 
involved in this self-inhibition may be released non synaptically from dendrites of 
the same or neighbouring cells making the calculation of this effect more complex 
than that for post-synaptic potentials. 
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